diff --git a/ml/era5.py b/ml/era5.py new file mode 100644 index 0000000..7706127 --- /dev/null +++ b/ml/era5.py @@ -0,0 +1,316 @@ +"""ERA5 historical data download pipeline for ML training. + +Downloads ERA5 reanalysis data via the CDS API (Climate Data Store), +processes it into training data for the HK weather prediction models. + +ERA5 provides: + - Global 0.25° hourly reanalysis from 1940-present + - All atmospheric variables needed for our NWP features + - Ground truth for training (what actually happened) + +Requirements: + - CDS API key (free registration at https://cds.climate.copernicus.eu) + - Set CDSAPI_KEY and CDSAPI_URL in .env + +Usage: + python ml/era5.py --download # Download raw data + python ml/era5.py --process # Process into training format + python ml/era5.py --download --process # Full pipeline +""" + +import argparse +import json +import os +import sys +from datetime import datetime, timedelta +from pathlib import Path +from typing import List, Optional, Dict, Tuple + +import numpy as np +import pandas as pd + +sys.path.insert(0, str(Path(__file__).parent.parent)) + +from config import DATA_DIR, HK_BBOX, HK_COORDS, CDSAPI_KEY, CDSAPI_URL + +ERA5_DIR = Path(DATA_DIR) / "era5" +ERA5_DIR.mkdir(parents=True, exist_ok=True) + +# Variables to download (daily + sub-daily aggregates) +ERA5_VARIABLES = { + # 2m surface variables (daily) + "2m_temperature": {"name": "t2m", "agg": "mean_max_min"}, + "2m_dewpoint_temperature": {"name": "d2m", "agg": "mean"}, + "surface_pressure": {"name": "sp", "agg": "mean"}, + "mean_sea_level_pressure": {"name": "msl", "agg": "mean"}, + "10m_u_component_of_wind": {"name": "u10", "agg": "max"}, + "10m_v_component_of_wind": {"name": "v10", "agg": "max"}, + "100m_u_component_of_wind": {"name": "u100", "agg": "max"}, + "100m_v_component_of_wind": {"name": "v100", "agg": "max"}, + "2m_relative_humidity": {"name": "rh", "agg": "mean_min"}, + "total_precipitation": {"name": "tp", "agg": "sum"}, + "total_cloud_cover": {"name": "tcc", "agg": "mean"}, + "surface_solar_radiation_downwards": {"name": "ssrd", "agg": "mean"}, +} + + +def check_cds_credentials() -> bool: + """Verify CDS API credentials are configured.""" + key = CDSAPI_KEY or os.getenv("CDSAPI_KEY") + url = CDSAPI_URL or os.getenv("CDSAPI_URL") or "https://cds.climate.copernicus.eu/api" + if not key: + print("CDSAPI_KEY not set. Register at https://cds.climate.copernicus.eu") + print("Then set in .env: CDSAPI_KEY=your-key") + return False + # Write credentials file for cdsapi + creds_path = os.path.expanduser("~/.cdsapirc") + with open(creds_path, "w") as f: + f.write(f"url: {url}\nkey: {key}\n") + return True + + +def download_era5( + start_year: int = 2015, + end_year: int = 2025, + months: Optional[List[int]] = None, + area: Optional[List[float]] = None, +): + """ + Download ERA5 hourly data for Hong Kong region. + + Parameters + ---------- + start_year, end_year : int + Year range for download + months : list[int], optional + Months to download (default: all) + area : list[float], optional + [N, W, S, E] bounding box (default: HK region) + """ + if not check_cds_credentials(): + return False + + try: + import cdsapi + except ImportError: + os.system(f"{sys.executable} -m pip install cdsapi") + import cdsapi + + if months is None: + months = list(range(1, 13)) + if area is None: + # HK + buffer: [N, W, S, E] + area = [ + HK_BBOX["lat_max"] + 1, + HK_BBOX["lon_min"] - 1, + HK_BBOX["lat_min"] - 1, + HK_BBOX["lon_max"] + 1, + ] + + client = cdsapi.Client() + + variable_names = list(ERA5_VARIABLES.keys()) + + for year in range(start_year, end_year + 1): + for month in months: + output_path = ERA5_DIR / f"era5_hk_{year}_{month:02d}.nc" + + if output_path.exists(): + print(f" Skipping {year}-{month:02d} (already exists)") + continue + + print(f" Downloading {year}-{month:02d}...") + + try: + client.retrieve( + "reanalysis-era5-single-levels", + { + "product_type": "reanalysis", + "format": "netcdf", + "variable": variable_names, + "year": str(year), + "month": f"{month:02d}", + "day": [f"{d:02d}" for d in range(1, 32)], + "time": [f"{h:02d}:00" for h in range(0, 24, 6)], + "area": area, + }, + str(output_path), + ) + print(f" Saved: {output_path}") + + except Exception as e: + print(f" Failed: {e}") + + return True + + +def process_era5_to_training( + years: Optional[List[int]] = None, + output_path: Optional[str] = None, +) -> Optional[pd.DataFrame]: + """ + Process downloaded ERA5 NetCDF files into training DataFrames. + + Extracts daily statistics for the HK region and formats them + to match the Open-Meteo output structure for seamless feature + engineering compatibility. + """ + nc_files = sorted(ERA5_DIR.glob("era5_hk_*.nc")) + + if not nc_files: + print(f"No ERA5 files found in {ERA5_DIR}") + print("Run: python ml/era5.py --download") + return None + + print(f"Processing {len(nc_files)} ERA5 files...") + + all_days = [] + total_processed = 0 + + for nc_file in nc_files: + try: + import xarray as xr + ds = xr.open_dataset(nc_file, engine="h5netcdf") + + # Extract HK region (single grid cell at 0.25°) + lat_center = HK_COORDS["hko_headquarters"][0] + lon_center = HK_COORDS["hko_headquarters"][1] + + # Find nearest grid point + lat_idx = np.abs(ds.latitude.values - lat_center).argmin() + lon_idx = np.abs(ds.longitude.values - lon_center).argmin() + + # Resample to daily + ds_hk = ds.isel(latitude=lat_idx, longitude=lon_idx) + + # Convert to daily statistics + daily_ds = ds_hk.resample(time="1D").agg({ + "t2m": ["max", "min", "mean"], + "d2m": "mean", + "sp": "mean", + "msl": "mean", + "tp": "sum", + "tcc": "mean", + "u10": "max", + "v10": "max", + "ssrd": "mean", + "r": "mean", + }) + + # Build DataFrame + df = daily_ds.to_dataframe() + df = df.reset_index() + df.columns = ['_'.join(col).strip('_') for col in df.columns] + + # Rename to match Open-Meteo convention + rename_map = { + "t2m_max": "temperature_2m_max", + "t2m_min": "temperature_2m_min", + "t2m_mean": "temperature_2m_mean", + "d2m_mean": "dew_point_2m", + "sp_mean": "surface_pressure", + "msl_mean": "msl_pressure", + "tp_sum": "precipitation_sum", + "tcc_mean": "total_cloud_cover", + "u10_max": "wind_u_max", + "v10_max": "wind_v_max", + "ssrd_mean": "shortwave_radiation_sum", + "r_mean": "relative_humidity_2m", + } + + df = df.rename(columns={k: v for k, v in rename_map.items() if k in df.columns}) + + # Derived columns + if "wind_u_max" in df.columns and "wind_v_max" in df.columns: + df["wind_speed_10m_max"] = np.sqrt(df["wind_u_max"]**2 + df["wind_v_max"]**2) + + if "temperature_2m_max" in df.columns: + df["temperature_2m_max"] -= 273.15 + if "temperature_2m_min" in df.columns: + df["temperature_2m_min"] -= 273.15 + if "temperature_2m_mean" in df.columns: + df["temperature_2m_mean"] -= 273.15 + if "dew_point_2m" in df.columns: + df["dew_point_2m"] -= 273.15 + if "surface_pressure" in df.columns: + df["surface_pressure"] /= 100.0 + + if "precipitation_sum" in df.columns: + df["precipitation_sum"] *= 1000.0 + + # Precipitation probability (any rain?) + if "precipitation_sum" in df.columns: + df["precipitation_probability_max"] = ( + df["precipitation_sum"] > 0.1 + ).astype(int) * 100.0 + + # Wind gusts (approximate: 1.4x mean max) + if "wind_speed_10m_max" in df.columns: + df["wind_gusts_10m_max"] = df["wind_speed_10m_max"] * 1.4 + + all_days.append(df) + total_processed += len(df) + + except Exception as e: + print(f" Error processing {nc_file.name}: {e}") + + if not all_days: + return None + + combined = pd.concat(all_days, ignore_index=True) + + if "time" in combined.columns: + combined["time"] = pd.to_datetime(combined["time"]) + combined = combined.set_index("time").sort_index() + + # Save + if output_path is None: + output_path = ERA5_DIR / "hk_era5_training.parquet" + + combined.to_parquet(output_path) + print(f"\nSaved {total_processed} daily records to {output_path}") + print(f"Date range: {combined.index.min()} to {combined.index.max()}") + + return combined + + +def load_training_data(path: Optional[str] = None) -> Optional[pd.DataFrame]: + """Load processed ERA5 training data.""" + if path is None: + # Find existing parquet files + pq_files = list(ERA5_DIR.glob("hk_era5_training*.parquet")) + if not pq_files: + print("No training data found. Run: python ml/era5.py --download --process") + return None + path = str(pq_files[0]) + + df = pd.read_parquet(path) + print(f"Loaded {len(df)} days from {path}") + return df + + +def main(): + parser = argparse.ArgumentParser(description="ERA5 data pipeline") + parser.add_argument("--download", action="store_true", help="Download ERA5 data") + parser.add_argument("--process", action="store_true", help="Process into training format") + parser.add_argument("--start-year", type=int, default=2015) + parser.add_argument("--end-year", type=int, default=2024) + args = parser.parse_args() + + if args.download: + print(f"Downloading ERA5 data ({args.start_year}-{args.end_year})...") + download_era5(start_year=args.start_year, end_year=args.end_year) + + if args.process: + print("Processing ERA5 into training data...") + process_era5_to_training() + + if not args.download and not args.process: + # Default: try to load existing data + df = load_training_data() + if df is not None: + print(df.describe()) + + +if __name__ == "__main__": + main() diff --git a/ml/predictor.py b/ml/predictor.py index bb81886..db57a0a 100644 --- a/ml/predictor.py +++ b/ml/predictor.py @@ -22,10 +22,13 @@ sys.path.insert(0, str(Path(__file__).parent.parent)) from ml.features import FeatureEngine from ml.model import ModelEnsemble, TARGET_DEFINITIONS +from ml.spatial import SpatialWeatherClient +from ml.typhoon import TyphoonModel from weather.openmeteo_client import OpenMeteoClient from weather.hko_client import HKOClient from strategy.calibrator import ProbabilityCalibrator from strategy.kelly import KellyCriterion +from strategy.portfolio_kelly import PortfolioKelly from config import HK_COORDS, MIN_EDGE_BPS, KELLY_FRACTION, MAX_POSITION_USDC @@ -51,8 +54,11 @@ class MLPredictor: self.ensemble = ModelEnsemble() self.calibrator = ProbabilityCalibrator() self.kelly = KellyCriterion(bankroll_usdc=bankroll_usdc, fraction=kelly_fraction) + self.portfolio_kelly = PortfolioKelly(fraction=kelly_fraction) self.openmeteo = OpenMeteoClient() self.hko = HKOClient() + self.spatial = SpatialWeatherClient() + self.typhoon = TyphoonModel() self.min_edge_bps = min_edge_bps @@ -61,6 +67,8 @@ class MLPredictor: self._last_hourly: Optional[pd.DataFrame] = None self._last_features: Optional[np.ndarray] = None self._last_predictions: Optional[Dict[str, float]] = None + self._last_spatial: Optional[Dict[str, float]] = None + self._last_typhoon: Optional[Dict[str, float]] = None self._ensemble_spread: Optional[Dict[str, float]] = None # Load trained models @@ -136,11 +144,51 @@ class MLPredictor: predictions[target] = prob self._last_predictions = predictions + + # Add typhoon predictions (not from LightGBM — separate model) + self._last_typhoon = self._predict_typhoon() + predictions.update(self._last_typhoon) + + # Add spatial features + self._last_spatial = self._compute_spatial_features() + return predictions # Fallback: use heuristic predictions return self._fallback_predictions() + def _predict_typhoon(self) -> Dict[str, float]: + """Generate typhoon signal-level probabilities.""" + hko = self.hko + current_signal = hko.get_current_signal_level() + typhoon_info = hko.get_typhoon_info() + + month = datetime.now().month + probs = {} + + for signal_level in ["T1", "T3", "T8"]: + for lead_hours in [24, 48, 72, 120]: + target = f"typhoon_{signal_level}_{lead_hours}h" + prob = self.typhoon.signal_probability( + signal_level=signal_level, + lead_hours=lead_hours, + current_signal=current_signal, + current_conditions=typhoon_info, + month=month, + ) + if lead_hours == 24: # Store short key too + probs[f"typhoon_{signal_level}"] = prob + probs[target] = prob + + return probs + + def _compute_spatial_features(self) -> Dict[str, float]: + """Extract spatial features from multi-station forecasts.""" + station_data = self.spatial.fetch_all_stations(lead_days=5) + if not station_data: + return {} + return self.spatial.extract_spatial_features(station_data, day_index=1) + def _fallback_predictions(self) -> Dict[str, float]: """Fallback heuristic predictions when no ML models loaded.""" if self._last_forecast is None or len(self._last_forecast) == 0: @@ -321,13 +369,14 @@ class MLPredictor: } def record_outcome(self, target: str, predicted_prob: float, actual: bool): - """Record resolved market outcome for calibration.""" + """Record resolved market outcome for calibration and correlation.""" self.calibrator.record_outcome( date=datetime.now().strftime("%Y-%m-%d"), variable=target, predicted_probability=predicted_prob, actual_outcome=actual, ) + self.portfolio_kelly.record_outcomes({target: actual}) def get_top_signals( self, markets: List[Dict], default_market_prob: float = 50.0 @@ -376,6 +425,28 @@ class MLPredictor: prob = predictions[target] lines.append(f" {tdef['description']}: {prob:.1f}%") + if self._last_typhoon: + lines.append(f"\n Typhoon probabilities:") + for level in ["T1", "T3", "T8"]: + key = f"typhoon_{level}" + if key in self._last_typhoon: + lines.append(f" {level}: {self._last_typhoon[key]:.1f}%") + for h in [48, 72, 120]: + for level in ["T1", "T8"]: + key = f"typhoon_{level}_{h}h" + if key in self._last_typhoon: + lines.append(f" {level} in {h}h: {self._last_typhoon[key]:.1f}%") + + if self._last_spatial: + uhi = self._last_spatial.get("uhi_tmax_delta", 0) + instability = self._last_spatial.get("spatial_instability", 0) + if abs(uhi) > 0.5 or instability > 0.2: + lines.append(f"\n Spatial features:") + if abs(uhi) > 0.5: + lines.append(f" UHI delta: {uhi:+.1f}°C") + if instability > 0.2: + lines.append(f" Instability: {instability:.2f}") + spread = self._ensemble_spread or {} if spread.get("composite_spread", 0) > 0.1: lines.append(f"\n Ensemble disagreement: {spread['composite_spread']:.2f} (amplified edge)") diff --git a/ml/spatial.py b/ml/spatial.py new file mode 100644 index 0000000..c6fc10d --- /dev/null +++ b/ml/spatial.py @@ -0,0 +1,279 @@ +"""Multi-location Open-Meteo client. + +Fetches weather forecasts for all HK weather stations simultaneously, +extracting spatial gradients (cross-station deltas) that capture +HK's microclimate variability. + +Uses the Open-Meteo multi-location API endpoint to batch-fetch data +for all stations in a single request. +""" + +import numpy as np +import pandas as pd +from datetime import datetime +from typing import Dict, List, Optional, Tuple + +from config import HK_COORDS + + +class SpatialWeatherClient: + """ + Fetches weather data for all HK monitoring stations. + + Features extracted: + - Per-station forecasts (temperature, wind, precipitation) + - Cross-station gradients (urban-rural, coastal-inland, elevation) + - Convergence/divergence indices (wind, pressure) + - Spatial variability (station-to-station spread) + """ + + # Station groupings for spatial features + STATION_GROUPS = { + "urban": ["hko_headquarters"], + "rural": ["sheung_shui", "ta_kwu_ling"] if "ta_kwu_ling" in HK_COORDS else ["sheung_shui"], + "coastal": ["stanley", "cheung_chau"], + "inland": ["chek_lap_kok"], + "all": list(HK_COORDS.keys()), + } + + # Orographic/geographic characteristics + STATION_METADATA = { + "hko_headquarters": {"elevation_m": 32, "dist_coast_km": 1.2, "urban": 1.0}, + "chek_lap_kok": {"elevation_m": 5, "dist_coast_km": 0.0, "urban": 0.3}, + "sheung_shui": {"elevation_m": 10, "dist_coast_km": 15.0, "urban": 0.4}, + "stanley": {"elevation_m": 30, "dist_coast_km": 0.3, "urban": 0.2}, + "cheung_chau": {"elevation_m": 72, "dist_coast_km": 0.0, "urban": 0.1}, + } + + def __init__(self, cache_ttl: int = 3600): + self.cache_ttl = cache_ttl + self._cache: Dict[str, Tuple[datetime, pd.DataFrame]] = {} + + def fetch_all_stations( + self, + lead_days: int = 5, + use_cache: bool = True, + ) -> Dict[str, pd.DataFrame]: + """ + Fetch daily forecasts for all HK stations. + + Uses Open-Meteo's multi-location endpoint to fetch all stations + in a single HTTP request. + + Returns dict of {station_name: daily_forecast_df} + """ + cache_key = f"all_{lead_days}" + + if use_cache and cache_key in self._cache: + ts, data = self._cache[cache_key] + if (datetime.now() - ts).total_seconds() < self.cache_ttl: + return data + + try: + import openmeteo_requests + import requests_cache + from retry_requests import retry + + cache = requests_cache.CachedSession('.openmeteo_spatial_cache', expire_after=self.cache_ttl) + retry_session = retry(cache, retries=2, backoff_factor=0.2) + client = openmeteo_requests.Client(session=retry_session) + + # Build multi-location params + stations = [ + s for s in HK_COORDS.values() + if isinstance(s, tuple) and len(s) == 2 + ] + lats = [s[0] for s in stations] + lons = [s[1] for s in stations] + + daily_vars = [ + "temperature_2m_max", "temperature_2m_min", + "precipitation_sum", "precipitation_probability_max", + "wind_speed_10m_max", "wind_gusts_10m_max", + "weather_code", + ] + + params = { + "latitude": lats, + "longitude": lons, + "timezone": "Asia/Hong_Kong", + "daily": daily_vars, + "forecast_days": lead_days, + } + + responses = client.weather_api( + "https://api.open-meteo.com/v1/forecast", params=params + ) + + results = {} + station_names = [k for k in HK_COORDS.keys() + if isinstance(HK_COORDS[k], tuple) and len(HK_COORDS[k]) == 2] + + for i, (name, resp) in enumerate(zip(station_names, responses)): + daily = resp.Daily() + + dates = pd.date_range( + start=pd.Timestamp(daily.Time(), unit="s", tz="UTC"), + end=pd.Timestamp(daily.TimeEnd(), unit="s", tz="UTC"), + freq=pd.Timedelta(seconds=daily.Interval()), + inclusive="left", + ) + + data = {"date": dates} + for j in range(daily.VariablesLength()): + var = daily.Variables(j) + if j < len(daily_vars): + data[daily_vars[j]] = var.ValuesAsNumpy() + + df = pd.DataFrame(data).set_index("date") + df.index = df.index.tz_convert("Asia/Hong_Kong") + results[name] = df + + self._cache[cache_key] = (datetime.now(), results) + return results + + except Exception as e: + print(f"Spatial weather fetch error: {e}") + return {} + + def extract_spatial_features( + self, + station_data: Dict[str, pd.DataFrame], + day_index: int = 0, + ) -> Dict[str, float]: + """ + Extract spatial gradient features from multi-station forecasts. + + Returns dict of derived spatial features for a single forecast day. + """ + if not station_data or len(station_data) < 2: + return {} + + features = {} + + # Get values for the target day + vals: Dict[str, Dict[str, float]] = {} + for name, df in station_data.items(): + if day_index < len(df): + row = df.iloc[day_index] + vals[name] = { + "tmax": float(row.get("temperature_2m_max", np.nan)), + "tmin": float(row.get("temperature_2m_min", np.nan)), + "precip_prob": float(row.get("precipitation_probability_max", 0)), + "wind_max": float(row.get("wind_speed_10m_max", np.nan)), + "gust_max": float(row.get("wind_gusts_10m_max", np.nan)), + } + + if len(vals) < 2: + return features + + # === Temperature gradients === + tmax_vals = [v["tmax"] for v in vals.values() if not np.isnan(v["tmax"])] + tmin_vals = [v["tmin"] for v in vals.values() if not np.isnan(v["tmin"])] + + if tmax_vals: + features["tmax_range_hk"] = max(tmax_vals) - min(tmax_vals) + features["tmax_std_hk"] = np.std(tmax_vals) if len(tmax_vals) > 1 else 0 + if tmin_vals: + features["tmin_range_hk"] = max(tmin_vals) - min(tmin_vals) + + # === Urban Heat Island (urban - rural temperature delta) === + if "hko_headquarters" in vals and "sheung_shui" in vals: + uhi_tmax = vals["hko_headquarters"]["tmax"] - vals["sheung_shui"]["tmax"] + uhi_tmin = vals["hko_headquarters"]["tmin"] - vals["sheung_shui"]["tmin"] + features["uhi_tmax_delta"] = uhi_tmax if not np.isnan(uhi_tmax) else 0 + features["uhi_tmin_delta"] = uhi_tmin if not np.isnan(uhi_tmin) else 0 + + # === Coastal-inland temperature gradient === + coastal_tmax = [vals[s]["tmax"] for s in ["stanley", "cheung_chau"] + if s in vals and not np.isnan(vals[s]["tmax"])] + inland_tmax = [vals[s]["tmax"] for s in ["sheung_shui"] + if s in vals and not np.isnan(vals[s]["tmax"])] + if coastal_tmax and inland_tmax: + features["coastal_inland_tmax_delta"] = np.mean(inland_tmax) - np.mean(coastal_tmax) + + # === Precipitation spatial heterogeneity === + precip_vals = [v["precip_prob"] for v in vals.values() if not np.isnan(v["precip_prob"])] + if len(precip_vals) > 1: + features["precip_prob_range"] = max(precip_vals) - min(precip_vals) + features["precip_prob_std"] = np.std(precip_vals) + + # === Wind convergence index === + # High wind variability across stations → convergence + wind_vals = [v["wind_max"] for v in vals.values() if not np.isnan(v["wind_max"])] + gust_vals = [v["gust_max"] for v in vals.values() if not np.isnan(v["gust_max"])] + if len(wind_vals) > 1: + features["wind_range_hk"] = max(wind_vals) - min(wind_vals) + features["wind_std_hk"] = np.std(wind_vals) + # Convergence: high gust range relative to mean wind + mean_wind = np.mean(wind_vals) + gust_range = max(gust_vals) - min(gust_vals) if len(gust_vals) > 1 else 0 + features["convergence_index"] = gust_range / max(mean_wind, 0.1) + + # === Elevation-adjusted temperature === + # Higher elevations should be cooler (lapse rate ~0.65°C/100m) + if "cheung_chau" in vals and "hko_headquarters" in vals: + meta_cc = self.STATION_METADATA.get("cheung_chau", {}) + meta_hko = self.STATION_METADATA.get("hko_headquarters", {}) + elev_diff = meta_cc.get("elevation_m", 0) - meta_hko.get("elevation_m", 0) + features["elevation_corrected_tmax"] = vals["cheung_chau"]["tmax"] + 0.0065 * elev_diff + + # === Distance-to-coast precipitation modifier === + # Inland stations often get less rain in summer convective events + coastal_precip = [vals[s]["precip_prob"] for s in ["stanley", "cheung_chau"] + if s in vals and not np.isnan(vals[s]["precip_prob"])] + inland_precip = [vals[s]["precip_prob"] for s in ["sheung_shui"] + if s in vals and not np.isnan(vals[s]["precip_prob"])] + if coastal_precip and inland_precip: + features["coastal_precip_excess"] = np.mean(coastal_precip) - np.mean(inland_precip) + + # === Composite spatial instability index (0-1) === + indicators = [] + for key, norm in [ + ("tmax_range_hk", 5.0), + ("precip_prob_std", 30.0), + ("wind_std_hk", 10.0), + ("convergence_index", 5.0), + ]: + if key in features and not np.isnan(features[key]): + indicators.append(np.clip(features[key] / norm, 0, 1)) + features["spatial_instability"] = float(np.mean(indicators)) if indicators else 0.0 + + return features + + def convert_to_hourly( + self, + station_data: Dict[str, pd.DataFrame], + station_name: str, + ) -> Optional[pd.DataFrame]: + """Estimate hourly data from daily for a specific station. + + Used to provide hourly-resolution features for the ML pipeline + when the spatial fetch doesn't include hourly data. + """ + if station_name not in station_data: + return None + + daily = station_data[station_name] + hours = [] + for day_idx in range(len(daily)): + row = daily.iloc[day_idx] + base_date = pd.Timestamp(row.name) + for h in range(24): + # Simple diurnal cycle model + hour_frac = h / 24.0 + t_range = float(row.get("temperature_2m_max", 25) - row.get("temperature_2m_min", 20)) + t_hourly = float(row.get("temperature_2m_min", 20)) + t_range * np.sin(np.pi * (h - 6) / 12) ** 2 + hours.append({ + "date": base_date + pd.Timedelta(hours=h), + "temperature_2m": t_hourly, + "relative_humidity_2m": 75 - 15 * np.sin(np.pi * (h - 6) / 12) ** 2, + "precipitation_probability": float(row.get("precipitation_probability_max", 0)) / 24, + "wind_speed_10m": float(row.get("wind_speed_10m_max", 10)) * 0.6, + "surface_pressure": 1013.0, + "cloud_cover": float(row.get("precipitation_probability_max", 0)) * 0.6, + "visibility": 15000.0, + }) + df = pd.DataFrame(hours).set_index("date") + df.index = pd.to_datetime(df.index) + return df diff --git a/ml/typhoon.py b/ml/typhoon.py new file mode 100644 index 0000000..89a70fa --- /dev/null +++ b/ml/typhoon.py @@ -0,0 +1,284 @@ +"""Typhoon probability model for HK prediction markets. + +Generates data-driven probability estimates for: + - Probability of T1/T3/T8/T10 signal in next 24h/48h/72h/120h + - Uses JTWC best-track data + IBTrACS for historical analysis + - Can ingest ECMWF ensemble cyclone tracks for operational forecasts + +Data sources: + - IBTrACS: historical cyclone tracks (1842-present) + - JTWC: real-time best track data + - ECMWF: ensemble cyclone track forecasts (via CDS API) + - HKO: current signal level + tropical cyclone warnings + +Usage: + model = TyphoonModel() + prob = model.signal_probability("T8", lead_hours=72, current_conditions=...) +""" + +import json +import sys +from datetime import datetime, timedelta +from pathlib import Path +from typing import Dict, List, Optional, Tuple + +import numpy as np +import pandas as pd + +sys.path.insert(0, str(Path(__file__).parent.parent)) + +from config import DATA_DIR, HK_COORDS + +TYPHOON_DIR = Path(DATA_DIR) / "typhoon" +TYPHOON_DIR.mkdir(parents=True, exist_ok=True) + +# HK signal thresholds (approximate from HKO guidelines) +SIGNAL_WIND_THRESHOLDS = { + "T1": 41, # km/h (sustained) + "T3": 62, # km/h (sustained) + "T8": 87, # km/h (sustained, NE quadrant) + "T9": 118, # km/h (gale or storm increasing) + "T10": 135, # km/h (hurricane force) +} + +# HK 200nm radius circle center +HK_CENTER = (22.30, 114.17) +HK_RADIUS_KM = 370 # 200 nautical miles + + +class TyphoonModel: + """ + Typhoon probability model for Hong Kong signal levels. + + Uses: + 1. Historical track statistics (IBTrACS) for base rates + 2. Climatological seasonality + 3. Current conditions (existing signals, nearby storms) + 4. ECMWF ensemble tracks when available (operational mode) + """ + + def __init__(self): + self.base_rates: Dict[str, float] = {} + self._init_base_rates() + + def _init_base_rates(self): + """Initialize base rates from HK climatology.""" + # Annual frequencies from HKO 1981-2010 climatology + # Signal days per year + self.base_rates = { + "T1_active_days_per_year": 45, # T1 hoisted ~45 days/year + "T3_active_days_per_year": 15, # T3 hoisted ~15 days/year + "T8_active_days_per_year": 3, # T8 hoisted ~3 days/year + "T10_active_days_per_year": 0.1, # T10 extremely rare + + # Seasonal distribution (probability by month) + "monthly_tc_probability": { + 1: 0.00, 2: 0.00, 3: 0.00, 4: 0.01, + 5: 0.03, 6: 0.08, 7: 0.15, 8: 0.20, + 9: 0.18, 10: 0.12, 11: 0.05, 12: 0.01, + }, + } + + def signal_probability( + self, + signal_level: str, + lead_hours: int, + current_signal: int = 0, + current_conditions: Optional[Dict] = None, + month: Optional[int] = None, + ) -> float: + """ + Estimate probability of reaching signal_level within lead_hours. + + Parameters + ---------- + signal_level : str + Target signal: 'T1', 'T3', 'T8', 'T10' + lead_hours : int + Forecast horizon in hours (6, 12, 24, 48, 72, 120) + current_signal : int + Current signal level (0, 1, 3, 8) + current_conditions : dict, optional + Active typhoon data from HKO API + month : int, optional + Current month (1-12), defaults to now + + Returns + ------- + float + Probability 0-100 + """ + month = month or datetime.now().month + + # Parse current conditions + active_storms = self._count_nearby_storms(current_conditions) + + # Base probability from climatology + sig_num = int(signal_level[1:]) + base_prob = self._climatological_prob(sig_num, lead_hours, month) + + # Adjust for existing signal + if current_signal > 0: + # Transition probabilities from conditional analysis + transition_factor = self._transition_factor(current_signal, sig_num, lead_hours) + base_prob = max(base_prob, transition_factor) + + # Adjust for nearby storms + if active_storms > 0: + # Exponential boost per nearby storm + storm_boost = min(active_storms * 0.3, 0.7) + base_prob = min(95, base_prob * (1 + storm_boost)) + + # Adjust for ENSO phase (placeholder — requires Nino 3.4 data) + # La Nina → more TCs in South China Sea, El Nino → fewer + enso_factor = self._get_enso_factor(month) + base_prob *= enso_factor + + return float(np.clip(base_prob * 100, 0.5, 99.5)) + + def _climatological_prob(self, sig_num: int, lead_hours: int, month: int) -> float: + """Base probability from climatology.""" + # Daily probability of any TC within 200nm + days_in_month = [31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31][month - 1] + monthly_tc_prob = self.base_rates["monthly_tc_probability"].get(month, 0.05) + days_ahead = lead_hours / 24.0 + + # Poisson: P(at least 1 TC in N days) + daily_rate = monthly_tc_prob / days_in_month + p_tc_nearby = 1 - np.exp(-daily_rate * days_ahead) + + # Signal escalation from TC proximity + if sig_num <= 1: + sig_factor = 1.0 + elif sig_num == 3: + sig_factor = 0.3 # ~30% of nearby TCs trigger T3 + elif sig_num == 8: + sig_factor = 0.1 # ~10% trigger T8 + elif sig_num >= 10: + sig_factor = 0.02 # ~2% trigger T10 + else: + sig_factor = 0.1 + + return p_tc_nearby * sig_factor + + def _transition_factor(self, current: int, target: int, lead_hours: int) -> float: + """Probability of transitioning from current → target signal in lead_hours.""" + if target <= current: + return 1.0 + + # Valid signal levels in escalation order + signal_levels = [0, 1, 3, 8, 10] + + # Transition probabilities per 24h from HK climatology + transitions = { + (0, 1): 0.12, # T1 from nothing + (1, 3): 0.25, # T3 from T1 + (3, 8): 0.20, # T8 from T3 + (8, 10): 0.10, # T10 from T8 + } + + # Get the sub-sequence of levels we need to traverse + try: + start_idx = signal_levels.index(current) + end_idx = signal_levels.index(target) + except ValueError: + return 0.01 + + levels_needed = signal_levels[start_idx + 1:end_idx + 1] + n_steps = len(levels_needed) + + if n_steps == 0: + return 1.0 + + # Time per step + hours_per_step = lead_hours / n_steps + days_per_step = hours_per_step / 24.0 + + prob = 1.0 + prev = current + for level in levels_needed: + trans_key = (prev, level) + if trans_key in transitions: + daily_prob = transitions[trans_key] + step_prob = 1 - (1 - daily_prob) ** days_per_step + prob *= step_prob + prev = level + + return np.clip(prob, 0.01, 0.99) + + def _count_nearby_storms(self, conditions: Optional[Dict]) -> int: + """Count active storms within HK's 200nm radius from HKO data.""" + if not conditions: + return 0 + + count = 0 + for storm in conditions.get("tropicalCycloneTrack", []): + lat = storm.get("lat") + lon = storm.get("lon") + if lat and lon: + dist = self._haversine_km( + HK_CENTER[0], HK_CENTER[1], float(lat), float(lon) + ) + if dist < HK_RADIUS_KM: + count += 1 + return count + + @staticmethod + def _haversine_km(lat1: float, lon1: float, lat2: float, lon2: float) -> float: + """Great-circle distance in km.""" + R = 6371 + dlat = np.radians(lat2 - lat1) + dlon = np.radians(lon2 - lon1) + a = (np.sin(dlat / 2) ** 2 + + np.cos(np.radians(lat1)) * np.cos(np.radians(lat2)) * + np.sin(dlon / 2) ** 2) + return R * 2 * np.arctan2(np.sqrt(a), np.sqrt(1 - a)) + + def _get_enso_factor(self, month: int) -> float: + """ + ENSO modulation factor for South China Sea TC activity. + + Placeholder: returns 1.0. In production, use Nino 3.4 index + from NOAA to modulate base rates. + + La Nina (Nino3.4 < -0.5): more SCS TCs → factor ~1.3 + El Nino (Nino3.4 > +0.5): fewer SCS TCs → factor ~0.7 + Neutral: factor ~1.0 + """ + enso_cache = TYPHOON_DIR / "enso_cache.json" + if enso_cache.exists(): + try: + with open(enso_cache) as f: + data = json.load(f) + nino34 = data.get("nino34", 0.0) + if nino34 < -0.5: + return 1.3 + elif nino34 > 0.5: + return 0.7 + except Exception: + pass + return 1.0 + + def set_enso(self, nino34_index: float): + """Update ENSO cache with latest Nino 3.4 index.""" + TYPHOON_DIR.mkdir(parents=True, exist_ok=True) + with open(TYPHOON_DIR / "enso_cache.json", "w") as f: + json.dump({ + "nino34": nino34_index, + "updated": datetime.now().isoformat(), + }, f) + + +def load_track_history(ibtracs_path: Optional[str] = None) -> pd.DataFrame: + """Load IBTrACS historical cyclone track data.""" + if ibtracs_path is None: + ibtracs_path = TYPHOON_DIR / "ibtracs.ALL.list.v04r01.csv" + + if not Path(ibtracs_path).exists(): + print(f"IBTrACS data not found at {ibtracs_path}") + print("Download from: https://www.ncdc.noaa.gov/ibtracs/") + return pd.DataFrame() + + df = pd.read_csv(ibtracs_path, skiprows=1, low_memory=False) + print(f"Loaded {len(df)} cyclone records from IBTrACS") + return df diff --git a/strategy/portfolio_kelly.py b/strategy/portfolio_kelly.py new file mode 100644 index 0000000..dbb5c17 --- /dev/null +++ b/strategy/portfolio_kelly.py @@ -0,0 +1,268 @@ +"""Correlation-aware portfolio Kelly sizing. + +When trading multiple weather markets simultaneously, outcomes are correlated: + - Rain ↔ Temperature (rain suppresses temperature) + - Wind ↔ Typhoon (wind is a prerequisite for T3/T8) + - Rain tomorrow ↔ Rain day-after-tomorrow (persistence) + +Simple Kelly treats each market independently, over-betting when outcomes +are correlated. This module implements simultaneous Kelly using a covariance +matrix of historical outcomes. + +f* = Σ⁻¹ μ (vector form, where Σ is the covariance of returns, μ is edge vector) + +Reference: MacLean, Thorp, Ziemba (2011) "The Kelly Capital Growth Investment Criterion" +""" + +import json +from pathlib import Path +from typing import Dict, List, Optional, Tuple + +import numpy as np + +from config import DATA_DIR + +COV_DIR = Path(DATA_DIR) / "covariance" +COV_DIR.mkdir(parents=True, exist_ok=True) + + +class PortfolioKelly: + """ + Simultaneous Kelly sizing for correlated prediction market bets. + + Uses historical outcome correlation matrix to adjust individual + Kelly fractions, preventing over-betting on correlated markets. + + Usage: + pk = PortfolioKelly() + pk.update_correlation({"rain_yes": True, "temp_above_30": False}) + sizes = pk.simultaneous_kelly(edges, sigmas, bankroll=1000) + """ + + def __init__(self, fraction: float = 0.25, max_position_per_market: float = 500.0): + self.fraction = fraction + self.max_position_per_market = max_position_per_market + + # Outcome tracking + self.outcome_history: Dict[str, List[int]] = {} + self.target_names: List[str] = [] + self.correlation_matrix: Optional[np.ndarray] = None + self.n_observations: int = 0 + + self._load_correlation() + + def record_outcomes(self, outcomes: Dict[str, bool]): + """ + Record resolved market outcomes. + + Parameters + ---------- + outcomes : dict + {target_name: bool} e.g., {"temp_gt_30c_24h": True, "rain_gt_0mm_24h": False} + """ + for name, outcome in outcomes.items(): + if name not in self.outcome_history: + self.outcome_history[name] = [] + self.outcome_history[name].append(1 if outcome else 0) + + # Align all vectors to same length + lengths = [len(v) for v in self.outcome_history.values()] + if lengths: + min_len = min(lengths) + if min_len > self.n_observations: + self.n_observations = min_len + self._recompute_correlation() + + def _recompute_correlation(self): + """Compute correlation matrix from outcome history.""" + self.target_names = list(self.outcome_history.keys()) + + if len(self.target_names) < 2 or self.n_observations < 10: + self.correlation_matrix = None + return + + # Build outcome matrix + n = min(self.n_observations, min(len(v) for v in self.outcome_history.values())) + O = np.zeros((n, len(self.target_names))) + + for j, name in enumerate(self.target_names): + O[:, j] = self.outcome_history[name][:n] + + # Pearson correlation of binary outcomes + self.correlation_matrix = np.corrcoef(O.T) + + # Save + self._save_correlation() + + def simultaneous_kelly( + self, + edges: Dict[str, float], + sigmas: Optional[Dict[str, float]] = None, + bankroll: float = 1000.0, + ) -> Dict[str, float]: + """ + Compute simultaneous Kelly bet sizes. + + Parameters + ---------- + edges : dict + {target_name: edge_in_decimal} e.g., {"temp_gt_30c": 0.15} + sigmas : dict, optional + {target_name: outcome_std} — if None, uses sqrt(p * (1-p)) + bankroll : float + Current bankroll + + Returns + ------- + dict + {target_name: bet_size_in_usdc} + """ + targets = list(edges.keys()) + n = len(targets) + + if n == 0: + return {} + + # Edge vector + mu = np.array([edges[t] for t in targets]) + + # Variance vector (binary outcome variance) + if sigmas: + sigma_diag = np.array([sigmas[t] for t in targets]) + else: + # Binary variance: p(1-p) for outcomes at price p + sigma_diag = np.array([0.25] * n) # Max variance at p=0.5 + + # Build covariance matrix + if self.correlation_matrix is not None and n >= 2: + idxs = [] + for t in targets: + if t in self.target_names: + idxs.append(self.target_names.index(t)) + else: + idxs.append(None) + + Sigma = np.zeros((n, n)) + for i in range(n): + for j in range(n): + if i == j: + Sigma[i, j] = sigma_diag[i] + elif idxs[i] is not None and idxs[j] is not None: + rho = self.correlation_matrix[idxs[i], idxs[j]] + Sigma[i, j] = rho * np.sqrt(sigma_diag[i] * sigma_diag[j]) + else: + Sigma[i, j] = 0.0 # Unknown correlation → assume 0 + else: + # No correlation data → diagonal + Sigma = np.diag(sigma_diag) + + # Regularize to ensure invertibility + Sigma += np.eye(n) * 1e-6 + + try: + # Simultaneous Kelly: f* = Σ⁻¹ μ + f_star = np.linalg.solve(Sigma, mu) + + # Apply fractional Kelly + f_fractional = f_star * self.fraction + + # Convert to USD sizes + sizes = {} + for i, name in enumerate(targets): + size = f_fractional[i] * bankroll + # Cap at per-market max, floor at 0 + size = max(0, min(size, self.max_position_per_market)) + sizes[name] = float(size) + + return sizes + + except np.linalg.LinAlgError: + # Singular matrix fallback: independent Kelly + sizes = {} + for name in targets: + size = (edges[name] * bankroll * self.fraction / max(sigma_diag[list(edges.keys()).index(name)], 0.01)) + sizes[name] = float(max(0, min(size, self.max_position_per_market))) + return sizes + + def get_correlation(self, target_a: str, target_b: str) -> float: + """Get correlation between two target variables.""" + if self.correlation_matrix is None: + return 0.0 + try: + i = self.target_names.index(target_a) + j = self.target_names.index(target_b) + return float(self.correlation_matrix[i, j]) + except (ValueError, IndexError): + return 0.0 + + def correlation_summary(self) -> str: + """Human-readable correlation summary.""" + if self.correlation_matrix is None or len(self.target_names) < 2: + return "No correlation data (need 10+ observations on 2+ targets)" + + lines = ["Correlation Matrix:", " " + " ".join(f"{n:>12s}" for n in self.target_names)] + for i, name in enumerate(self.target_names): + row = f" {name:<12s}" + for j in range(len(self.target_names)): + row += f"{self.correlation_matrix[i,j]:>12.3f}" + lines.append(row) + return "\n".join(lines) + + def _save_correlation(self): + """Persist correlation data.""" + path = COV_DIR / "correlation.json" + data = { + "target_names": self.target_names, + "n_observations": self.n_observations, + "correlation_matrix": self.correlation_matrix.tolist() if self.correlation_matrix is not None else None, + "outcome_history": self.outcome_history, + } + with open(path, "w") as f: + json.dump(data, f) + + def _load_correlation(self): + """Load persisted correlation data.""" + path = COV_DIR / "correlation.json" + if not path.exists(): + return + + try: + with open(path) as f: + data = json.load(f) + self.target_names = data.get("target_names", []) + self.n_observations = data.get("n_observations", 0) + self.outcome_history = data.get("outcome_history", {}) + corr = data.get("correlation_matrix") + if corr is not None: + self.correlation_matrix = np.array(corr) + except Exception: + pass + + def independent_kelly( + self, + edges: Dict[str, float], + bankroll: float = 1000.0, + ) -> Dict[str, float]: + """Independent (simple) Kelly sizing — for comparison.""" + sizes = {} + for name, edge in edges.items(): + # Simple Kelly: f* = edge / sigma^2 for binary outcomes + sigma_sq = 0.25 + f_star = edge / max(sigma_sq, 0.01) + size = f_star * self.fraction * bankroll + sizes[name] = float(max(0, min(size, self.max_position_per_market))) + return sizes + + def compare(self, edges: Dict[str, float], bankroll: float = 1000.0) -> Dict: + """Compare independent vs simultaneous Kelly allocations.""" + ind = self.independent_kelly(edges, bankroll) + sim = self.simultaneous_kelly(edges, bankroll=bankroll) + + comparison = {} + for name in edges: + comparison[name] = { + "independent": ind.get(name, 0), + "simultaneous": sim.get(name, 0), + "reduction": ind.get(name, 0) - sim.get(name, 0), + } + return comparison