from pathlib import Path
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from joblib import Memory
from IPython.display import display
from cycler import cycler
from scipy.stats import norm
from quantfinlab.common.cache import cache_key
from quantfinlab.dataio.realtime import read_statscan_vintages, read_alfred
from quantfinlab.dataio.macro import read_macro_forecasts, read_boc_market
from quantfinlab.macro import realtime, transforms, bridge, inflation, midas, dfm, bvar, policy, evaluation
from quantfinlab.ml.combination import forecast_weights, combine_forecasts
from quantfinlab.ml.probabilistic import draw_summary, gaussian_crps, gaussian_nll, interval_coverage
from quantfinlab.fixed_income.overnight import reference_window, compound_overnight
from quantfinlab.plotting import macro as macro_plots
palette = ["#069AF3", "#FE420F", "#00008B", "#008080", "#CC79A7", "#9614fa",
"#DC143C", "#7BC8F6", "#0072B2", "#04D8B2", "#800080", "#FF8072"]
plt.rcParams["axes.prop_cycle"] = cycler(color=palette)
plt.rcParams.update({"figure.figsize": (6, 3), "figure.dpi": 200, "savefig.dpi": 300,
"axes.grid": True, "grid.alpha": .20, "axes.spines.top": False,
"axes.spines.right": False, "axes.titlesize": 12, "axes.labelsize": 12,
"xtick.labelsize": 9, "ytick.labelsize": 9, "legend.fontsize": 7})
data = Path("../data")
source = data / "canada_statscan_realtime.parquet"
definitions = {
"monthly_gdp": ("36-10-0491", "Canada | Seasonally adjusted at annual rates | Chained (2017) dollars | All industries [T001]"),
"manufacturing_gdp": ("36-10-0491", "Canada | Seasonally adjusted at annual rates | Chained (2017) dollars | Manufacturing [31-33]"),
"employment": ("14-10-0331", "Canada | Employment for all employees | Industrial aggregate excluding unclassified businesses [11-91N]"),
"earnings": ("14-10-0331", "Canada | Average weekly earnings including overtime for all employees | Industrial aggregate excluding unclassified businesses [11-91N]"),
"manufacturing_sales": ("16-10-0014", "Canada | Sales of goods manufactured (shipments) | Total, durable and non-durable goods"),
"manufacturing_orders": ("16-10-0014", "Canada | New orders | Total, durable and non-durable goods"),
"manufacturing_inventories": ("16-10-0014", "Canada | Inventories | Total, durable and non-durable goods"),
"wholesale_volume": ("20-10-0005", "Canada | Wholesale sales Chained Fisher volume index (scaled to equal 100 in 2012) | Wholesale trade [41]"),
"retail_volume": ("20-10-0082", "Canada | Retail trade [44-45] | Retail sales at 2012 constant prices (Unchained Laspeyres Index)"),
"goods_exports": ("12-10-0165", "Canada | Export | Balance of payments | Seasonally adjusted | Total of all merchandise"),
"goods_imports": ("12-10-0165", "Canada | Import | Balance of payments | Seasonally adjusted | Total of all merchandise"),
"core_trim": ("18-10-0259", "Canada | Measure of core inflation based on a trimmed mean approach, CPI-trim (year-over-year percent change)"),
"core_median": ("18-10-0259", "Canada | Measure of core inflation based on a weighted median approach, CPI-median (year-over-year percent change)"),
"core_common": ("18-10-0259", "Canada | Measure of core inflation based on a factor model, CPI-common (year-over-year percent change)")}
components = {"consumption": "Household final consumption expenditure", "residential": "Residential structures",
"business": "Non-residential structures, machinery and equipment", "ip": "Intellectual property products",
"exports": "Exports of goods and services", "imports": "Less: imports of goods and services",
"government": "General governments final consumption expenditure", "real_gdp": "Gross domestic product at market prices"}
for name, estimate in components.items():
definitions[name] = ("36-10-0431", f"Canada | Chained (2017) dollars | Seasonally adjusted at annual rates | {estimate}")
definitions[f"nominal_{name}"] = ("36-10-0431", f"Canada | Current prices | Seasonally adjusted at annual rates | {estimate}")
selected = {name: {"table_id": table, "series_title": title} for name, (table, title) in definitions.items()}
vintages = read_statscan_vintages(source, selected, start="1995-01-01")
cpi_definitions = {name: {"table_id": "18-10-0004", "series_title": f"Canada | {title}"}
for name, title in {"headline_cpi": "All-items", "core_cpi": "All-items excluding food and energy",
"food_cpi": "Food", "energy_cpi": "Energy"}.items()}
cpi = read_statscan_vintages(data / "canada_statscan_current.parquet", cpi_definitions, start="1995-01-01")
cpi_calendar = realtime.first_releases(vintages[vintages["series_id"].eq("core_trim")], max_delay=100).set_index("observation_date")["available_at"]
cpi["available_at"] = cpi["observation_date"].map(cpi_calendar)
cpi["timing"] = np.where(cpi["available_at"].notna(), "CPI release matched to core vintage", "next-month-end availability bound")
cpi["available_at"] = cpi["available_at"].fillna(cpi["observation_date"] + pd.offsets.MonthEnd(2))
vintages = pd.concat([vintages, cpi], ignore_index=True)
us = read_alfred(data / "alfred_realtime.parquet", series=["INDPRO", "RSAFS", "IPMAN", "PAYEMS", "UNRATE", "ICSA", "CCSA", "PCEC96"], start="1995-01-01")
us["series_id"] = "us_" + us["series_id"]
bos = read_macro_forecasts(data / "canada_boc_bos.csv")
bos = bos[bos["series_id"].isin(["FUTURESALES", "INVESTMACHEQUIP", "EMPLOY", "CREDIT", "INPUTS", "OUTPUTS"])
& bos["release_date"].notna()].rename(columns={"observation_quarter": "observation_date", "release_date": "available_at"})
bos["series_id"] = "bos_" + bos["series_id"]
boc = read_boc_market(data / "canada_boc_market.parquet")
boc_daily = boc.pivot(index="date", columns="series_id", values="value").sort_index()
controls = boc[boc["series_id"].isin(["target_overnight", "usd_cad", "goc_2y", "goc_5y", "goc_10y", "corra"])].rename(columns={"date": "observation_date"})
controls["available_at"] = controls["observation_date"]
controls["series_id"] = controls["series_id"].replace({"target_overnight": "policy"})
information = pd.concat([vintages, us, bos, controls], ignore_index=True)
quarterly_names = list(components) + [f"nominal_{name}" for name in components]
monthly_names = [name for name in definitions if name not in quarterly_names]
monthly_names += list(cpi_definitions) + [f"us_{name}" for name in ["INDPRO", "RSAFS", "IPMAN", "PAYEMS", "UNRATE", "ICSA", "CCSA", "PCEC96"]]
monthly_names += list(bos["series_id"].unique()) + ["policy", "usd_cad", "goc_2y", "goc_5y", "goc_10y"]
information = information[information["series_id"].isin(monthly_names)][
["series_id", "observation_date", "available_at", "valid_to", "value"]].reset_index(drop=True)
frequencies = {name: "Q" for name in bos["series_id"].unique()}
codes = {name: 5 for name in monthly_names}
codes.update({name: 1 for name in ["core_trim", "core_median", "core_common", *frequencies]})
codes.update({name: 2 for name in ["policy", "goc_2y", "goc_5y", "goc_10y", "us_UNRATE"]})
codes.update({name: 6 for name in cpi_definitions})
key = cache_key([source, data / "canada_statscan_current.parquet", data / "alfred_realtime.parquet",
data / "canada_boc_market.parquet", data / "canada_boc_bos.csv"],
{"calculation": "canada-nowcasting-library-1", "series": selected, "codes": codes, "seed": 23})
cache = Path("../workspace/library_repeat/cache") / key
memory = Memory(cache, verbose=0)
snapshot = realtime.monthly_snapshot
release_growth = memory.cache(realtime.release_growth)
fit_factors = memory.cache(dfm.fit_dfm)
targets = {"real_gdp": ("Q", 1, 100, "percent"), "headline_cpi": ("M", 12, 100, "percent"),
"core_trim": ("M", 1, 1, "level"), "core_median": ("M", 1, 1, "level"),
"employment": ("M", 12, 100, "percent"), "earnings": ("M", 12, 100, "percent")}
truth_parts = []
for name, (frequency, periods, scale, method) in targets.items():
source_values = vintages[vintages["series_id"].eq(name)]
values = release_growth(source_values, periods=periods, frequency=frequency, scale=scale,
method=method, max_delay=180 if frequency == "Q" else 100,
annualization=4 if frequency == "Q" else 1)
latest_levels = realtime.vintage_asof(source_values, source_values["available_at"].max(), wide=True)[name]
latest_growth = transforms.growth_rate(latest_levels, periods=periods, scale=scale, method=method,
annualization=4 if frequency == "Q" else 1)
values["latest"] = values["observation_date"].map(latest_growth)
values["target"] = name
values["exact_release"] = values["observation_date"].isin(cpi_calendar.index) if name == "headline_cpi" else True
values["unit"] = "Annualized q/q (%)" if frequency == "Q" else "Year-over-year (%)"
truth_parts.append(values)
truth = pd.concat(truth_parts, ignore_index=True).dropna(subset=["first"]).sort_values(["target", "observation_date"])
days = {name: [1, 7, 15, 30, 45, 60] if name == "real_gdp" else [1, 5, 10, 20] for name in targets}
grid = evaluation.forecast_origins(truth[truth["exact_release"]], days, start="2015-01-01")
score_start = pd.Timestamp("2019-01-01")
gdp_grid = grid[grid["target"].eq("real_gdp")]
component_truth = []
for name in components:
values = release_growth(vintages[vintages["series_id"].eq(name)], periods=1,
frequency="Q", scale=100, method="percent", annualization=4, max_delay=180)
component_truth.append(values.assign(component=name))
component_truth = pd.concat(component_truth).pivot(index="observation_date", columns="component", values="first")
blocks = {"consumption": ["monthly_gdp", "wholesale_volume", "retail_volume", "us_RSAFS"],
"residential": ["monthly_gdp", "policy", "goc_5y"],
"business": ["manufacturing_sales", "manufacturing_orders", "us_IPMAN"],
"ip": ["monthly_gdp", "us_INDPRO"], "exports": ["goods_exports", "us_INDPRO", "usd_cad"],
"imports": ["goods_imports", "manufacturing_sales", "usd_cad"],
"government": ["employment", "monthly_gdp"],
"common_activity": ["monthly_gdp", "manufacturing_inventories", "us_INDPRO", "us_PAYEMS"]}
gdp_rows = []
for row in gdp_grid.itertuples():
known = snapshot(information, row.evaluation_date, frequencies=frequencies, series=monthly_names, start="1997-01-01")
record = row._asdict()
for name in dict.fromkeys(name for block in blocks.values() for name in block):
signal = bridge.quarterly_signal(known[name].dropna(), row.observation_date.to_period("Q"),
difference=name in ["policy", "goc_5y"])
record[name], record[f"{name}_known"] = signal["signal"], signal["observed_months"]
nominal = realtime.vintage_asof(vintages, row.evaluation_date, series=[f"nominal_{name}" for name in components], wide=True)
nominal = nominal[nominal.index < row.observation_date].dropna().iloc[-1]
for component in blocks:
if component != "common_activity":
record[f"w_{component}"] = nominal[f"nominal_{component}"] / nominal["nominal_real_gdp"] * (-1 if component == "imports" else 1)
gdp_rows.append(record)
gdp_design = pd.DataFrame(gdp_rows).merge(component_truth, on="observation_date", suffixes=("", "_component"))
gdp_design["common_activity"] = gdp_design["actual"] - sum(
gdp_design[f"w_{name}"] * gdp_design[name] for name in blocks if name != "common_activity")
gdp_design["w_common_activity"] = 1.0
bridge_rows, contribution_rows = [], []
for row in gdp_design.itertuples(index=False):
current = pd.Series(row._asdict())
train = gdp_design[gdp_design["horizon"].eq(row.horizon) & gdp_design["release_date"].lt(row.evaluation_date)]
contributions, variances = {}, []
for component, names in blocks.items():
prediction, sigma, _ = bridge.bridge_nowcast(train, current, [*names, *[f"{name}_known" for name in names]], component, alpha=5, minimum=20)
contributions[component] = current[f"w_{component}"] * prediction
variances.append((current[f"w_{component}"] * sigma) ** 2)
bridge_rows.append({"target": "real_gdp", "observation_date": row.observation_date, "release_date": row.release_date,
"evaluation_date": row.evaluation_date, "horizon": row.horizon, "actual": row.actual,
"model": "Component bridge", "mean": sum(contributions.values()), "sigma": np.sqrt(np.nansum(variances))})
contribution_rows.append({"observation_date": row.observation_date, "horizon": row.horizon, **contributions})
forecasts = [pd.DataFrame(bridge_rows)]
contributions = pd.DataFrame(contribution_rows)
factor_names = ["global", "activity", "labor", "inflation", "financial"]
groups = {}
for name in monthly_names:
group = "inflation" if name in [*cpi_definitions, "core_trim", "core_median", "core_common", "bos_INPUTS", "bos_OUTPUTS"] else "activity"
if name in ["employment", "earnings", "us_PAYEMS", "us_UNRATE", "us_ICSA", "us_CCSA", "bos_EMPLOY"]:
group = "labor"
if name in ["policy", "usd_cad", "goc_2y", "goc_5y", "goc_10y", "bos_CREDIT"]:
group = "financial"
groups[name] = ["global", group]
groups["real_gdp"] = ["global", "activity"]
orders = {"global": 2, "activity": 1, "labor": 1, "inflation": 1, "financial": 1}
anchors = {"global": "monthly_gdp", "activity": "monthly_gdp", "labor": "employment", "inflation": "core_trim", "financial": "goc_5y"}
factor_fits, states, systems, baseline_rows, factor_rows, var_rows = {}, {}, {}, [], [], []
for row in grid[grid["release_date"].ge(score_start)].itertuples(index=False):
date = row.evaluation_date
target_history = truth[truth["target"].eq(row.target) & truth["release_date"].lt(date)].sort_values("observation_date")
if target_history.empty:
continue
record = {"target": row.target, "observation_date": row.observation_date, "release_date": row.release_date,
"evaluation_date": date, "horizon": row.horizon, "actual": row.actual}
for model, prediction in evaluation.baseline_forecast(target_history["first"], quarterly=row.target == "real_gdp").items():
baseline_rows.append({**record, "model": model, "mean": prediction, "sigma": target_history["first"].tail(60).std(ddof=1)})
if len(target_history) < (20 if row.target == "real_gdp" else 36):
continue
regime = pd.Timestamp(2017 + 2 * ((date.year - 2017) // 2), 1, 1)
if regime not in factor_fits:
raw = snapshot(information, regime, frequencies=frequencies, series=monthly_names, start="1997-01-01")
transformed = transforms.fred_transform(raw, codes)
usable = transformed.columns[transformed.notna().sum().ge(48)].tolist()
quarterly = truth[truth["target"].eq("real_gdp") & truth["release_date"].le(regime)].set_index("observation_date")[["first"]].rename(columns={"first": "real_gdp"})
factor_fits[regime] = fit_factors(transformed[usable], quarterly,
factors={name: groups[name] for name in [*usable, "real_gdp"]}, orders=orders, anchors=anchors)
fit = factor_fits[regime]
if date not in states:
raw = snapshot(information, date, frequencies=frequencies, series=monthly_names, start="1997-01-01")
transformed = transforms.fred_transform(raw, codes)
quarterly = truth[truth["target"].eq("real_gdp") & truth["release_date"].le(date)].set_index("observation_date")[["first"]].rename(columns={"first": "real_gdp"})
states[date] = {"factors": dfm.filter_dfm(fit, transformed, quarterly)["factors"]}
factors = states[date]["factors"]
train = dfm.factor_history(factors, target_history, frequency=targets[row.target][0])
current = dfm.factor_point(factors, row.observation_date, frequency=targets[row.target][0])
current["last_release"] = target_history["first"].iloc[-1]
prediction = dfm.factor_forecast(train, current, [*factor_names, "last_release"], "first")
factor_rows.append({**record, "model": "Grouped DFM", "mean": prediction["mean"], "sigma": prediction["sigma"]})
if date not in systems:
monthly = truth[truth["release_date"].le(date) & truth["target"].ne("real_gdp")].pivot(index="observation_date", columns="target", values="first")
policy_monthly = boc_daily.loc[:date, "target_overnight"].resample("MS").mean().rename("policy")
monthly = factors[["activity"]].join(monthly).join(policy_monthly).dropna()
gaps = monthly.index.to_period("M").asi8
first = np.r_[0, np.flatnonzero(np.diff(gaps) != 1) + 1][-1]
monthly = monthly.iloc[first:]
if len(monthly) >= 40:
fitted = bvar.fit_bvar(monthly, persistent=["policy"])
paths = bvar.bvar_paths(fitted, steps=18, draws=300, rng=np.random.default_rng(23))
systems[date] = fitted, paths
if date in systems:
fitted, paths = systems[date]
if row.target == "real_gdp":
draws = bvar.quarterly_bridge_draws(fitted, paths, factors["activity"], target_history,
row.observation_date.to_period("Q"), rng=np.random.default_rng(23))
else:
step = max(1, (row.observation_date.to_period("M") - fitted["z"].index.max().to_period("M")).n)
draws = paths[:, min(step, paths.shape[1]) - 1, fitted["columns"].index(row.target)]
var_rows.append({**record, "model": "Minnesota BVAR", **draw_summary(draws)})
baselines = pd.DataFrame(baseline_rows).sort_values("release_date")
errors = baselines["mean"] - baselines["actual"]
baselines["sigma"] = errors.groupby([baselines["target"], baselines["horizon"], baselines["model"]]).transform(
lambda x: x.expanding(min_periods=12).std(ddof=1).shift(1))
forecasts.extend([baselines, pd.DataFrame(factor_rows), pd.DataFrame(var_rows)])
high_frequency = read_macro_forecasts(data / "macro_high_frequency.parquet")
high_frequency = high_frequency.pivot(index="date", columns="series_id", values="value").sort_index()
oil = high_frequency["DCOILBRENTEU"].dropna()
gas = high_frequency["GASREGW"].dropna()
fx = boc_daily["usd_cad"].dropna()
oil = oil * fx.reindex(oil.index, method="ffill")
gas = gas * fx.reindex(gas.index, method="ffill")
weekly_oil = 100 * np.log(oil.resample("W-FRI").last()).diff()
weekly_gas = 100 * np.log(gas.resample("W-FRI").last()).diff()
component_history = []
for name, label in [("core_cpi", "core"), ("food_cpi", "food"), ("energy_cpi", "energy"), ("headline_cpi", "headline")]:
values = release_growth(vintages[vintages["series_id"].eq(name)], periods=12,
scale=100, method="percent", max_delay=100)
component_history.append(values[["observation_date", "release_date", "first"]].rename(columns={"first": label}))
components_monthly = component_history[0]
for values in component_history[1:]:
components_monthly = components_monthly.merge(values, on=["observation_date", "release_date"])
energy_design = []
for row in components_monthly.itertuples():
for h in days["headline_cpi"]:
date = row.release_date - pd.offsets.BDay(h)
energy_design.append({"observation_date": row.observation_date, "horizon": h,
"gas_signal": inflation.market_growth(gas, date, row.observation_date, periods=12, scale=100),
"oil_signal": inflation.market_growth(oil, date, row.observation_date, periods=12, scale=100)})
energy_design = pd.DataFrame(energy_design)
component_rows, bridge_inflation_rows, lag_rows = [], [], []
for row in grid[grid["release_date"].ge(score_start)].itertuples(index=False):
record = {"target": row.target, "observation_date": row.observation_date, "release_date": row.release_date,
"evaluation_date": row.evaluation_date, "horizon": row.horizon, "actual": row.actual}
if row.target == "headline_cpi":
history = components_monthly[components_monthly["release_date"].lt(row.evaluation_date)].copy()
history["last_energy"] = history["energy"].shift(1)
train = history.merge(energy_design[energy_design["horizon"].eq(row.horizon)], on="observation_date")
current = pd.Series({"last_energy": history["energy"].iloc[-1],
"gas_signal": inflation.market_growth(gas, row.evaluation_date, row.observation_date, periods=12, scale=100),
"oil_signal": inflation.market_growth(oil, row.evaluation_date, row.observation_date, periods=12, scale=100)})
result = inflation.inflation_components(history, train, current, energy="energy", prior=(.70, .17, .13))
past = pd.DataFrame(component_rows)
if len(past):
past = past[past["horizon"].eq(row.horizon) & past["release_date"].lt(row.evaluation_date)]
sigma = (past["mean"] - past["actual"]).tail(60).std(ddof=1) if len(past) >= 12 else history["headline"].tail(36).std(ddof=1)
component_rows.append({**record, "model": "CPI components", "mean": result["mean"],
"sigma": sigma,
**result["components"].to_dict(),
**{f"{name}_contribution": result["components"][name] * result["weights"][name]
for name in ["core", "food", "energy"]}})
if row.target in ["core_trim", "core_median"]:
known = realtime.vintage_asof(vintages, row.evaluation_date, series=["headline_cpi", "core_cpi"], wide=True)
change = transforms.growth_rate(known, periods=12, method="percent")
history = truth[truth["target"].eq(row.target) & truth["release_date"].lt(row.evaluation_date)].copy()
train = history.merge(components_monthly[["observation_date", "core", "headline"]], on="observation_date")
train["core_lag"] = train["core"].shift(1)
train["headline_lag"] = train["headline"].shift(1)
train["last_release"] = train["first"].shift(1)
train["current_available"] = False
current = pd.Series({"core_lag": change["core_cpi"].dropna().iloc[-1],
"headline_lag": change["headline_cpi"].dropna().iloc[-1],
"last_release": history["first"].iloc[-1], "current_available": False})
fit = inflation.fit_inflation_bridge(train, current, features=["core_lag", "headline_lag", "last_release"], target="first")
prediction = inflation.inflation_forecast(fit, current)
bridge_inflation_rows.append({**record, "model": "Inflation bridge", **prediction})
if row.target in ["headline_cpi", "employment", "real_gdp"]:
record = record.copy()
if row.target == "headline_cpi":
signals = {"oil": midas.release_lags(weekly_oil, row.evaluation_date, 12),
"gas": midas.release_lags(weekly_gas, row.evaluation_date, 12)}
elif row.target == "employment":
claims = realtime.vintage_asof(us, row.evaluation_date, series=["us_ICSA"], wide=True)["us_ICSA"].dropna()
signals = {"claims": midas.release_lags(100 * np.log(claims).diff(), row.evaluation_date, 12)}
elif row.evaluation_date in states:
signals = {"activity": midas.release_lags(states[row.evaluation_date]["factors"]["activity"], row.evaluation_date, 6)}
else:
signals = {}
for name, values in signals.items():
record.update({f"{name}_{lag}": value for lag, value in enumerate(values)})
lag_rows.append(record)
lags = pd.DataFrame(lag_rows)
midas_rows = []
for target, names in {"headline_cpi": {"oil": 12, "gas": 12}, "employment": {"claims": 12}, "real_gdp": {"activity": 6}}.items():
design = lags[lags["target"].eq(target)].sort_values("release_date")
weights = {}
shape_cutoff = pd.Timestamp("2021-01-01")
for name, length in names.items():
columns = [f"{name}_{lag}" for lag in range(length)]
train = design[design["release_date"].lt(shape_cutoff)].dropna(subset=columns)
shape, _ = midas.fit_beta_shape(train[columns].to_numpy(), train["actual"].to_numpy())
weights[name] = midas.beta_weights(length, *shape)
for row in design[design["evaluation_date"].ge(shape_cutoff)].itertuples(index=False):
train = design[design["horizon"].eq(row.horizon) & design["release_date"].lt(row.evaluation_date)].copy()
train["last_release"] = train["actual"].shift(1)
columns = [f"{name}_{lag}" for name, length in names.items() for lag in range(length)]
train = train.dropna(subset=["actual", "last_release", *columns])
minimum = 12 if target == "real_gdp" else 36
if len(train) >= minimum and all(pd.notna(getattr(row, column)) for column in columns):
fitted = midas.fit_midas({name: train[[f"{name}_{lag}" for lag in range(length)]].to_numpy() for name, length in names.items()},
train["actual"], controls=train[["last_release"]], weights=weights, minimum=minimum)
current = {name: np.array([[getattr(row, f"{name}_{lag}") for lag in range(length)]]) for name, length in names.items()}
prediction = midas.midas_forecast(fitted, current, controls=pd.DataFrame({"last_release": [train["actual"].iloc[-1]]}))[0]
midas_rows.append({"target": target, "observation_date": row.observation_date, "release_date": row.release_date,
"evaluation_date": row.evaluation_date, "horizon": row.horizon, "actual": row.actual,
"model": "MIDAS", "mean": prediction, "sigma": fitted["regression"].sigma})
forecasts.extend([pd.DataFrame(component_rows), pd.DataFrame(bridge_inflation_rows), pd.DataFrame(midas_rows)])
forecasts = pd.concat(forecasts, ignore_index=True).dropna(subset=["mean"])
forecasts = forecasts[forecasts["release_date"].ge(score_start)].sort_values(["release_date", "target", "horizon", "model"])
averages, weight_rows = [], []
for (target, observation_date, date, h), candidates in forecasts.groupby(["target", "observation_date", "evaluation_date", "horizon"], sort=True):
candidates = candidates[candidates["model"].ne("Rolling mean")]
history = forecasts[forecasts["target"].eq(target) & forecasts["horizon"].eq(h)]
values = truth[truth["target"].eq(target) & truth["release_date"].lt(date)]["first"].tail(60)
scale = max(1.4826 * (values - values.median()).abs().median(), values.std(ddof=1) / 3, 1e-6)
weights = forecast_weights(history, date, scale=scale)
weights = weights.reindex(candidates["model"]).dropna()
if weights.empty:
fallback = candidates.loc[candidates["model"].isin(["Last release", "AR(1)"]), "model"]
weights = pd.Series(1 / len(fallback), index=fallback)
weights /= weights.sum()
candidates = candidates.set_index("model").loc[weights.index]
mean = candidates["mean"].clip(values.median() - 8 * scale, values.median() + 8 * scale)
sigma = candidates["sigma"].fillna(scale).clip(lower=.1 * scale)
result = combine_forecasts(mean, sigma, weights)
averages.append({"target": target, "observation_date": observation_date, "evaluation_date": date,
"release_date": candidates["release_date"].iloc[0], "horizon": h,
"actual": candidates["actual"].iloc[0], "model": "Adaptive average", **result,
"q10": result["mean"] + norm.ppf(.1) * result["sigma"],
"q90": result["mean"] + norm.ppf(.9) * result["sigma"]})
for name, weight in weights.items():
weight_rows.append({"target": target, "observation_date": observation_date,
"evaluation_date": date, "horizon": h, "model": name, "weight": weight})
forecasts = pd.concat([forecasts, pd.DataFrame(averages)], ignore_index=True)
weights = pd.DataFrame(weight_rows)
policy_rows, corra_rows = [], []
system_dates = pd.Series(sorted(systems)).groupby(pd.DatetimeIndex(sorted(systems)).to_period("M")).max()
for date in system_dates:
fitted, paths = systems[date]
monthly = fitted["z"] * fitted["scale"] + fitted["center"]
for h in [3, 6, 12]:
target_month = date.to_period("M") + h
target_date = target_month.to_timestamp()
actual = boc_daily.loc[boc_daily.index.to_period("M") == target_month, "target_overnight"].mean()
step = (target_month - monthly.index.max().to_period("M")).n
if 0 < step <= paths.shape[1]:
rule = policy.fit_policy_rule(monthly, features=["activity", "core_trim", "employment", "policy"], horizon=h)
prediction = policy.policy_forecast(rule, monthly.iloc[[-1]])
draws = paths[:, step - 1, fitted["columns"].index("policy")]
policy_rows.append({"origin": date, "target_date": target_date, "horizon": h,
"actual": actual, "BVAR": draws.mean(), "Policy rule": prediction["mean"][0],
"Random walk": monthly["policy"].iloc[-1]})
next_quarter = date.to_period("Q") + 1
start, end = reference_window(pd.date_range(next_quarter.start_time,
next_quarter.start_time + pd.offsets.MonthEnd(0), freq="W-WED")[2])
fixings = boc_daily["corra"].dropna() / 100
calendar = fixings.index[(fixings.index >= start) & (fixings.index < end)]
if end <= fixings.index.max() and len(calendar) and calendar[0] == start:
months = pd.date_range(monthly.index.max() + pd.offsets.MonthBegin(1), periods=paths.shape[1], freq="MS")
basis = ((boc_daily["corra"] - boc_daily["target_overnight"]).loc[:date].dropna().tail(60).median()) / 100
draws = policy.policy_window_draws(paths[:, :, fitted["columns"].index("policy")] / 100,
months, fixings, as_of=date, start=start, end=end,
calendar=calendar, basis=basis, day_count=365)
actual = compound_overnight(fixings.reindex(calendar).to_numpy(), calendar, end, day_count=365)
corra_rows.append({"origin": date, "start": start, "end": end, "actual": actual * 100,
"BVAR": draws.mean() * 100, "Random walk": fixings.loc[:date].iloc[-1] * 100})
policy_history = pd.DataFrame(policy_rows).dropna(subset=["actual"])
corra_history = pd.DataFrame(corra_rows).sort_values("origin").drop_duplicates("start", keep="last")
policy_scores = []
for h, group in policy_history.groupby("horizon"):
for model in ["BVAR", "Policy rule", "Random walk"]:
error = group[model] - group["actual"]
policy_scores.append({"comparison": f"Policy mean, {h}m", "model": model, "n": len(error),
"RMSE": np.sqrt(np.mean(error ** 2)), "MAE": error.abs().mean()})
for model in ["BVAR", "Random walk"]:
error = corra_history[model] - corra_history["actual"]
policy_scores.append({"comparison": "Three-month compounded CORRA", "model": model, "n": len(error),
"RMSE": np.sqrt(np.mean(error ** 2)), "MAE": error.abs().mean()})
policy_scores = pd.DataFrame(policy_scores)
mps = read_macro_forecasts(data / "canada_boc_mps.csv")
survey = mps[mps["question"].str.startswith("1.7") & mps["row_label"].eq("Median of responses")
& mps["column_label"].str.contains("End of")].copy()
survey["year"] = survey["column_label"].str.extract(r"End of (\d{4})").astype(int)
survey_rows = []
for row in survey.itertuples():
eligible = system_dates[system_dates.le(row.release_date)]
if eligible.empty:
continue
date = eligible.iloc[-1]
fitted, paths = systems[date]
month = pd.Period(f"{row.year}-12", freq="M")
step = (month - fitted["z"].index.max().to_period("M")).n
actual = truth[truth["target"].eq("headline_cpi") & truth["observation_date"].eq(month.to_timestamp())]
if 0 < step <= paths.shape[1] and not actual.empty:
survey_rows.append({"release_date": row.release_date, "origin": date, "target_date": month.to_timestamp(),
"actual": actual["first"].iloc[0], "MPS median": row.value_numeric,
"BVAR": paths[:, step - 1, fitted["columns"].index("headline_cpi")].mean()})
survey_comparison = pd.DataFrame(survey_rows)
survey_scores = []
for model in ["MPS median", "BVAR"]:
error = survey_comparison[model] - survey_comparison["actual"]
survey_scores.append({"comparison": "Year-end CPI; MPS public-release comparison", "model": model,
"n": len(error), "RMSE": np.sqrt(np.mean(error ** 2)), "MAE": error.abs().mean()})
policy_scores = pd.concat([policy_scores, pd.DataFrame(survey_scores)], ignore_index=True)
after_date = max(states)
before_date = max(date for date in states if date <= after_date - pd.Timedelta(days=30)
and date.year // 2 == after_date.year // 2)
regime = pd.Timestamp(2017 + 2 * ((after_date.year - 2017) // 2), 1, 1)
fit = factor_fits[regime]
impact_date = after_date.to_period("Q").end_time.to_period("M").to_timestamp()
news_states = []
for date in [before_date, after_date]:
raw = snapshot(information, date, frequencies=frequencies, series=monthly_names, start="1997-01-01")
quarterly = truth[truth["target"].eq("real_gdp") & truth["release_date"].le(date)].set_index("observation_date")[["first"]].rename(columns={"first": "real_gdp"})
news_states.append(dfm.filter_dfm(fit, transforms.fred_transform(raw, codes), quarterly, end=impact_date)["result"])
news = dfm.dfm_news(news_states[0], news_states[1], variable="real_gdp",
impact_date=impact_date, location=fit.quarterly_location["real_gdp"], scale=fit.quarterly_scale["real_gdp"])
news_summary = news["impacts"].set_index(["impact date", "impacted variable"])
release_rows = []
for target, model in [("headline_cpi", "CPI components"), ("employment", "MIDAS"), ("earnings", "AR(1)")]:
sample = forecasts[forecasts["target"].eq(target) & forecasts["model"].eq(model)
& forecasts["horizon"].eq(1)].sort_values("release_date").copy()
sample["news"] = evaluation.release_surprises(sample["actual"].reset_index(drop=True), sample["mean"].reset_index(drop=True), minimum=24, cap=5).to_numpy()
release_rows.append(sample[["observation_date", "release_date", "news"]].rename(columns={"news": target}))
employment_news = release_rows[1].merge(release_rows[2], on=["observation_date", "release_date"])
yield_history = boc_daily[["goc_2y", "goc_5y", "goc_10y"]].dropna() / 100
reaction = pd.concat([
evaluation.release_response(release_rows[0], yield_history, drivers=["headline_cpi"]).assign(release="CPI"),
evaluation.release_response(employment_news, yield_history, drivers=["employment", "earnings"]).assign(release="SEPH bundle")])
matched = []
for (target, h), group in forecasts.groupby(["target", "horizon"]):
wide = group.pivot(index="observation_date", columns="model", values="mean")
common = wide.dropna().index
matched.append(group[group["observation_date"].isin(common)])
matched = pd.concat(matched, ignore_index=True)
scores = evaluation.nowcast_scores(matched)
score_table = scores.pivot(index=["target", "horizon"], columns="model", values="RMSE")
coverage = truth.groupby("target").agg(first_observation=("observation_date", "min"),
last_observation=("observation_date", "max"), releases=("first", "size"), unit=("unit", "first"))
coverage["scored_first_release"] = matched.groupby("target")["release_date"].min()
coverage["scored_last_release"] = matched.groupby("target")["release_date"].max()
coverage["scored_observations"] = matched.groupby("target")["observation_date"].nunique()
latest = forecasts[forecasts["horizon"].eq(1)].sort_values("evaluation_date").groupby(["target", "model"]).tail(1)
latest = latest.pivot(index=["target", "observation_date", "evaluation_date"], columns="model", values="mean")
density_rows = []
for target, group in forecasts[forecasts["model"].eq("Minnesota BVAR")].groupby("target"):
group = group.dropna(subset=["actual", "mean", "sigma", "q10", "q90"])
density_rows.append({"target": target, "NLL": gaussian_nll(group["actual"], group["mean"], group["sigma"] ** 2),
"Gaussian CRPS": np.mean(gaussian_crps(group["actual"], group["mean"], group["sigma"])),
"80% coverage": interval_coverage(group["actual"], group["q10"], group["q90"]),
"interval width": (group["q90"] - group["q10"]).mean()})
density = pd.DataFrame(density_rows).set_index("target")
response_table = reaction[reaction["market"].eq("goc_2y")].copy()
response_table.index = response_table["release"] + " / " + response_table["driver"]
diagnostics = pd.concat({"Density": density, "2Y yield response": response_table[["beta", "se", "p_value", "events"]]})
display(coverage,
score_table.style.format(precision=3, na_rep="—").set_caption("RMSE on common observations within each target and horizon"),
latest.style.format(precision=3, na_rep="—").set_caption("Last historical evaluation; one business day before release"),
news_summary.round(4), policy_scores.set_index(["comparison", "model"]).round(3),
diagnostics.style.format(precision=4, na_rep="—").set_caption("BVAR densities: all available origins; yield response: close-to-close"))
loadings = dfm.factor_loadings(fit)
labels = {"monthly_gdp": "Monthly GDP", "manufacturing_gdp": "Manufacturing GDP",
"manufacturing_sales": "Manufacturing sales", "manufacturing_orders": "Manufacturing orders",
"core_trim": "CPI-trim", "core_median": "CPI-median", "core_common": "CPI-common",
"employment": "SEPH employment", "earnings": "SEPH earnings", "us_PAYEMS": "US payroll",
"us_UNRATE": "US unemployment", "us_ICSA": "US initial claims", "us_CCSA": "US continuing claims",
"us_RSAFS": "US retail sales", "us_INDPRO": "US industrial production", "us_IPMAN": "US manufacturing",
"goc_2y": "Canada 2Y yield", "goc_5y": "Canada 5Y yield", "goc_10y": "Canada 10Y yield",
"bos_INPUTS": "BOS: input prices", "bos_OUTPUTS": "BOS: selling prices",
"bos_EMPLOY": "BOS: future employment", "bos_INVESTMACHEQUIP": "BOS: investment",
"bos_FUTURESALES": "BOS: future sales", "bos_CREDIT": "BOS: credit conditions"}
loadings = loadings.rename(index=labels)
news_details = news["details"].copy()
news_details["updated variable"] = news_details["updated variable"].replace(labels)
gdp = forecasts[forecasts["target"].eq("real_gdp") & forecasts["horizon"].eq(7)
& forecasts["model"].isin(["Component bridge", "Grouped DFM", "Minnesota BVAR", "AR(1)"])]
weight_history = weights[weights["target"].eq("real_gdp") & weights["horizon"].eq(30)].pivot(
index="evaluation_date", columns="model", values="weight")
inflation_history = pd.DataFrame(component_rows)
inflation_history = inflation_history[inflation_history["horizon"].eq(1)].set_index("observation_date")[["actual", "mean", "core_contribution", "food_contribution", "energy_contribution"]]
inflation_history.columns = ["Headline CPI", "Component forecast", "Core contribution", "Food contribution", "Energy contribution"]
last_policy_date = max(systems)
fitted, paths = systems[last_policy_date]
future_months = pd.date_range(fitted["z"].index.max() + pd.offsets.MonthBegin(1), periods=paths.shape[1], freq="MS")
policy_path = paths[:, :, fitted["columns"].index("policy")]
policy_plot = pd.DataFrame({"mean": policy_path.mean(axis=0), "q10": np.quantile(policy_path, .1, axis=0),
"q90": np.quantile(policy_path, .9, axis=0)}, index=future_months)
policy_plot["actual"] = boc_daily["target_overnight"].resample("MS").mean().reindex(future_months)
policy_plot = policy_plot[policy_plot.index.to_period("M") > last_policy_date.to_period("M")].head(12)
fig, axes = plt.subplots(4, 2, figsize=(16, 20), constrained_layout=True)
macro_plots.plot_nowcast(gdp, ax=axes[0, 0], title="Canadian GDP: seven business days before release")
axes[0, 0].set_yscale("symlog", linthresh=3, linscale=.8)
axes[0, 0].set_yticks([-100, -40, -10, -3, 0, 3, 10, 40], labels=["−100", "−40", "−10", "−3", "0", "3", "10", "40"])
axes[0, 0].set_ylabel("Annualized q/q % (symmetric log)")
macro_plots.plot_nowcast_scores(scores[scores["target"].eq("real_gdp")], ax=axes[0, 1])
macro_plots.plot_revisions(truth[truth["target"].eq("real_gdp")], ax=axes[1, 0])
macro_plots.plot_factor_loadings(loadings.fillna(0), ax=axes[1, 1])
macro_plots.plot_inflation_components(inflation_history, ax=axes[2, 0], unit="Percentage points (y/y)")
axes[2, 0].set_title("Headline inflation and component contributions")
macro_plots.plot_forecast_weights(weight_history, ax=axes[2, 1])
macro_plots.plot_news_impacts(news_details, ax=axes[3, 0])
macro_plots.plot_policy_path(policy_plot, ax=axes[3, 1], title=f"BoC outlook at {last_policy_date:%Y-%m-%d}\nLast complete BVAR month: {fitted['z'].index.max():%Y-%m}")
fig.suptitle("Canada: real-time macro and monetary policy · library repeat", fontsize=16)
plt.show()