Time series econometrics: from the correlogram to the Kalman filter¶

This guide follows the order of a graduate time-series course (the chapter numbers are those of Neusser, Time Series Econometrics, Springer 2016) and names the sp call for each step. Every function here has a reference check against R or Stata on committed data; the table at the end says which.

import numpy as np
import pandas as pd
import statspai as sp

One series¶

Autocorrelation and the long-run variance (chapters 4 and 5)¶

sp.corrgram(df, "growth", lags=12)        # ACF, PACF, Ljung-Box Q
sp.ljungbox(df["growth"], lags=[4, 8])    # portmanteau test alone
sp.lrvar(df, "growth")                    # sum of all autocovariances

The variance of a sample mean of an autocorrelated series is the long-run variance divided by T, not the variance divided by T. sp.lrvar returns both, with the kernel and the bandwidth it used. The default is the Bartlett kernel with Andrews' data-driven bandwidth. kernel="qs", prewhite=1 is the choice of the sandwich package in R.

fit = sp.arima(df["growth"], order=(1, 0, 3))
fit.summary()
fit.forecast(8)
sp.arima(df["growth"], auto=True, stepwise=False)   # search over orders

Likelihoods of mixed ARMA models have several maxima. When two orders give nearly the same fit, prefer the smaller one, and do not read much into the coefficients of a model whose AR and MA roots nearly cancel.

The frequency domain and filters (chapter 6)¶

spec = sp.periodogram(df, "growth", method="smoothed", spans=[5, 5])
spec.table.head()          # freq, cycle_length, spectrum, lower, upper
spec.plot()
sp.periodogram(df, "growth", method="ar")           # AR spectral estimate
sp.cumulative_periodogram_test(df, "resid")         # Bartlett's white-noise test

The raw periodogram (method="raw") is not a consistent estimate of the spectral density at any single frequency. Smooth it, or fit the AR spectrum.

hp = sp.tsfilter(df, "lgdp", method="hp", smooth=1600)
bk = sp.tsfilter(df, "lgdp", method="bk", low=6, high=32, K=12)
hp.trend, hp.cycle
hp.gain(np.linspace(0, np.pi, 200))    # what the filter does to each frequency

method= also takes "cf" (Christiano-Fitzgerald), "bw" (Butterworth) and "hamilton" (the regression of y[t + h] on p lags). sp.tsfilter(df, "lgdp", one_sided=True) is the real-time HP trend, which uses data up to each date only and is the one to use when the cycle feeds a forecast or a policy rule. Three cautions. The HP filter applied to data that are not seasonally adjusted leaves the seasonal component in the cycle. The HP filter is two-sided, so the last few cycle values are revised when new data arrive. And any of these filters will produce a smooth "cycle" from a random walk; the gain function shows which frequencies it lets through, not whether the data have power there.

Trends, unit roots and breaks (chapter 7)¶

sp.unitroot(df, "lgdp", test="adf", trend="ct")     # also "pp", "kpss", "dfgls"
sp.zivot_andrews(df, "lgdp", model="both", lags=2)  # one break, date unknown
bn = sp.beveridge_nelson("lgdp", data=df)            # AR(p) by OLS
sp.beveridge_nelson("lgdp", data=df, order=(1, 1))  # ARMA by exact ML
bn.trend, bn.cycle, bn.long_run_multiplier

A break in the level or the slope of a trend looks like a unit root to a test that does not allow for it. sp.zivot_andrews picks the break date least favourable to the unit root and uses critical values that account for the search. A rejection supports "stationary around a trend that breaks once"; it does not date the break with any stated precision.

The Beveridge-Nelson trend is the long-horizon forecast of the series. Its cycle depends on the AR order of the differenced series, sometimes a lot; look at two or three orders before interpreting it.

Volatility (chapter 8)¶

sp.estat(ols_fit, "archlm", lags=5)                 # ARCH LM test
fit = sp.garch("ret", data=df, p=1, q=1, ar=1, dist="t")
fit.summary()
fit.forecast(10)           # conditional variance path
fit.value_at_risk(0.01)    # 1% quantile of the next observation

p counts lagged variances and q lagged squared innovations. ar= adds autoregressive terms to the mean and dist="t" estimates the degrees of freedom of a Student t innovation.

sp.garch("ret", data=df, model="gjr")       # threshold GARCH
sp.garch("ret", data=df, model="egarch")    # exponential GARCH

in_mean=True puts the conditional variance in the mean (a risk premium) and works with every model; in_mean="sd" uses the standard deviation and in_mean_lags=[0, 1] adds a lag.

Both asymmetric models let bad news move the variance by more than good news. In the threshold model gamma > 0 is that leverage effect; in EGARCH it is theta < 0. Stata's tarch coefficient has the opposite sign of gamma because Stata attaches it to positive shocks. A coefficient estimated at zero is on the boundary of the parameter space, where the usual standard errors do not apply; the fit warns, and the lower-order model is the one to report.

Several series¶

Cross-correlation (chapter 11)¶

cc = sp.xcorr("gdp", "sentiment", data=df, lags=8, prewhiten="ar")
cc.table      # lag, xcorr, band, outside
cc.haugh      # portmanteau test of no relation at any lead or lag

At lag h the table holds the correlation of gdp[t + h] with sentiment[t], so a large value at a positive lag says that sentiment moves first. Without prewhiten=, the autocorrelation of each series leaks into every lag and the table cannot be read.

VAR, structural VAR and their bands (chapters 12 to 15)¶

sp.varsoc(df, maxlag=8)                 # lag order
fit = sp.var(df, lags=2)
fit.granger_table(); fit.stability(); fit.lm_test(4)
fit.forecast(8)

out = sp.irf(fit, periods=20, ci="asymptotic")        # delta method
out = sp.irf(fit, periods=20, ci="bootstrap", reps=1000, seed=1)
out["irf"]["x -> y"], out["lower"]["x -> y"], out["upper"]["x -> y"]
fit.fevd(20, ci="asymptotic")        # variance shares with standard errors

Orthogonalised responses depend on the order of the variables. The bootstrap regenerates the series from resampled residuals and re-estimates the VAR. The default band is bias-corrected (boot="kilian"): the bootstrap replicates of a VAR are biased towards less persistence, and with a root near one in a short sample the plain percentile band covers far less often than it claims. boot="efron" is that plain band, the one Stata and R's vars report; it warns when the largest estimated root exceeds 0.9.

nan = np.nan
bq = sp.svar(fit, long_run=[[nan, 0], [nan, nan]])    # Blanchard-Quah
bq.irf(40, cumulative=True, ci="bootstrap", reps=500, seed=1)
ab = sp.svar(fit, A=A_pattern, B=B_pattern)           # short-run AB model
sg = sp.svar(fit, sign={"demand": {"y": "+", "p": "+"}})

Bands of an AB or long-run model repeat the identification on every bootstrap sample. Sign restrictions give a set of models, not one; their irf() already reports the spread of that set, which is not a confidence band.

An external instrument identifies one shock without restricting the others, and every point-identified model decomposes the data history:

px = sp.svar(fit, proxy=df["news"], proxy_target="i")   # proxy SVAR
px.proxy_diagnostics["robust_f"]          # proxy strength; warns below 10
px.irf(20, ci="bootstrap", reps=500, seed=1)   # moving block bootstrap of (u, proxy) pairs
px.fevd(20)                               # share of the full forecast-error variance
ab.historical_decomposition()             # actual = baseline + shock contributions

Local projections: IV, state-dependent and panel¶

sp.local_projections(df, outcome="y", shock="g", instruments=["news"],
                     lags=4, endog_order=["y", "g"], horizons=16)    # lpirfs::lp_lin_iv
sp.local_projections(df, outcome="y", shock="g", lags=4, state="gap",
                     gamma=3, hp_lambda=1600, horizons=16)          # lpirfs::lp_nl_iv
sp.local_projections(panel, outcome="y", shock="x", id="country",
                     time="year", lags=2, horizons=8)               # lpirfs::lp_lin_panel
sp.local_projections(df, outcome="y", shock="g", lags=2,
                     endog_order=["y"], common_sample=True,
                     vce="ols", dfk=False, horizons=7)              # Stata lpirf, exog(g)
sp.granger_causality(data=panel, caused="y", causing="x", lags=1,
                     id="country", time="year")                     # Dumitrescu-Hurlin

Any of lags, instruments, vce, trend, common_sample or dfk selects the explicit design: lags of the endog_order variables replace auto_lag, and the Newey-West lag at horizon h is h + 1 (the lpirfs rule) unless nw_lags is given. IV projections report the first-stage F of every horizon and warn below 10. State-dependent projections return the responses in the two regimes of the weight F and a per-horizon test of their equality. Panel projections cluster by unit by default; lpirfs reports homoskedastic errors unless asked, vce="ols" reproduces them.

Cointegration (chapter 16)¶

rank = sp.johansen(df, lags=1, trend="c")             # trace test
vec = sp.vec(df, lags=1, rank=2)
sp.johansen_lrtest(df, rank=2, beta_known=[1, 0, -1, 0])   # c - y stationary?
sp.johansen_lrtest(df, rank=2, beta=H)                     # beta = H phi
sp.johansen_lrtest(df, rank=2, loading=np.eye(4)[:, :3])   # weak exogeneity

lags counts lagged differences: Stata's vecrank, lags(2) and ca.jo(K = 2) are lags=1 here. The rank tests have non-standard distributions. Once the rank is fixed, hypotheses about the cointegrating vectors and the loadings are ordinary chi-squared tests, and they are where the economics usually sits.

State space models and the Kalman filter (chapter 17)¶

Write the model as X[t] = F X[t-1] + V[t], Y[t] = A + G X[t] + W[t] with Var(V) = Q and Var(W) = R.

res = sp.kalman_filter(y, F=F, G=G, Q=Q, R=R)     # known matrices
res.states("smoothed"); res.forecast(8)

def build(theta):                                  # unknown parameters
    return {"F": theta[0], "G": 1.0,
            "Q": np.exp(theta[1]), "R": np.exp(theta[2])}

fit = sp.statespace(y, build, start=[0.5, 0.0, 0.0])
fit.table; fit.filter.states("smoothed")

Missing values in y are allowed anywhere, which is what makes mixed frequencies easy: a quarterly state observed through an annual series is an observation equation that is missing three quarters out of four. Any matrix may vary over time (give it a leading time axis). For a regression with slowly moving coefficients sp.dlm is shorter. init="exact" starts non-stationary states from an exactly diffuse prior; the default for them is a large finite variance.

Regimes (chapter 18)¶

ms = sp.mswitch("growth", data=df, states=2, model="ar", ar=1)   # Hamilton
ms.summary()
ms.transition, ms.durations
ms.smoothed          # probability of each regime at each date
sp.mswitch("growth", data=df, states=2, switch_variance=True)   # dr model

model="ar" is Hamilton's model, in which the series adjusts gradually to a new regime mean. model="dr" lets the mean jump at once and treats lags of the series as ordinary regressors. The likelihood has several maxima, so several starting values are tried and the result says how many distinct maxima they reached.

test = sp.mswitch_lrtest("growth", data=df, model="ar", ar=1, reps=199, seed=1)
test.pvalue, test.critical_values, test.diagnostics

A likelihood-ratio test of one regime against two does not have a chi-squared distribution, because the transition probabilities are not identified under the null. sp.mswitch_lrtest simulates the null distribution from the fitted one-regime model instead. It takes minutes, not seconds; n_jobs= spreads the replicates over processes.

Coefficients that move (chapter 18)¶

tv = sp.tvp_var(df, variables=["dy", "dp", "r"], lags=2)   # Kalman, ML variances
tv.coefficients().head()      # date, equation, term, coef, se, band
tv.stability()                # largest root of the VAR frozen at each date
tv.irf(at=["1975Q1", "2005Q1"], periods=12)
sp.tvp_var(df, lags=2, method="forgetting", lam=0.99, kappa=0.96)

sv = sp.tvp_var_sv(df, variables=["dp", "dy", "r"], lags=2, seed=1)
sv.volatility()                       # standard deviation of each shock over time
sv.irf(at=["1975Q1", "2005Q1"], periods=12, shock_size="unit")
sv.diagnostics()

sp.tvp_var_sv also lets the shock variances drift, by Gibbs sampling. Use it when the question is whether responses changed or only the size of the shocks did: a model with constant variances has to explain a quieter period with different coefficients. Its default prior on coefficient drift is tight; check the sensitivity to k_Q.

Every coefficient follows a random walk. method="kalman" estimates how fast each one moves; a variance estimated at zero means that coefficient is constant. method="forgetting" discounts old data at a fixed rate and estimates nothing. The two can disagree about whether a response changed over time. Report both before claiming that it did.

What each function is checked against¶

Function	Reference
`sp.corrgram`, `sp.ljungbox`	Stata `corrgram`, R `Box.test`
`sp.lrvar`	R `sandwich`
`sp.arima`	R `forecast`, Stata `arima`
`sp.periodogram`	R `spec.pgram`, `spec.ar`; Stata `pergram`
`sp.cumulative_periodogram_test`	Stata `wntestb`
`sp.tsfilter`	Stata `tsfilter`, statsmodels
`sp.unitroot`	statsmodels, Stata
`sp.zivot_andrews`	R `urca::ur.za`, Stata `zandrews`
`sp.beveridge_nelson`	its definition (no package found)
`sp.garch`	Stata `arch`, R `rugarch`
`sp.xcorr`	R `ccf`, Stata `xcorr`
`sp.var`, `sp.irf`, `sp.svar`	Stata `var`, `irf`, `svar`; R `vars`
`sp.svar(proxy=)`	the 2SLS identity on the residuals; recovery of the true impact column
`SVARResult.historical_decomposition`	R `svars::hd`; the shock-free path from the initial values
`sp.local_projections` (IV, state, panel)	R `lpirfs` `lp_lin_iv`, `lp_nl_iv`, `lp_nl`, `lp_lin_panel`; Stata `lpirf`
`sp.granger_causality(id=, time=)`	R `plm::pgrangertest`, Stata `xtgcause`
`sp.johansen`, `sp.vec`, `sp.johansen_lrtest`	Stata `vecrank`, `vec`; R `urca`
`sp.kalman_filter`, `sp.statespace`	R `KFAS`, statsmodels
`sp.mswitch`	Stata `mswitch`, statsmodels
`sp.mswitch_lrtest`	rejection rates in simulations (no reference implementation)
`sp.tvp_var`	`sp.dlm` and R `KFAS` (Kalman); discounted least squares (forgetting)
`sp.tvp_var_sv`	exact checks of each Gibbs block and a joint-distribution test; R `bvarsv` as a screen

The bootstrap bands are random and have no cross-language reference; they are checked by coverage in simulations.