Dynamic Regression Models
2026-07-22
Screenshot
B-J modeling strategy starts with (p, d, q) diagnostics:
e_{t} = \phi_{1} e_{t-1} + \phi_{2}e_{t-2} + \cdots + \phi_{p} e_{t-p} + \nu_{t} + \theta_{1} \nu_{t-1} + \theta_{2}\nu_{t-2} + \cdots + \theta_{q} \nu_{t-q}
Determine Stationarity through a three-pronged strategy:
Diagnose the AR and the MA processes
Post-estimation diagnostics of the residuals from a filtering model — if we have the B-J model right, \hat{\epsilon} should be pure random noise.
Goal: White noise residuals – that is, residuals that contain no notable autocorrelation.
Measure: Portmanteau tests are encompassing tests for white noise in a model’s residuals. One of the most popular is the Ljung-Box Q-statistic.
If Y is white noise, then (because correlations are distributed normal, their squares are chi-square):
Q = n(n+2)\sum_{j=1}^{m} \frac{1}{n-j}\hat{\rho}_j^2 \;\rightarrow\; \chi^2_m
Note
Larger Q-statistics indicate that the errors contain some type of systematic process.
corrgram Output — An ExampleLAG AC PAC Q Prob>Q
1 0.9458 0.9706 92.164 0.0000
2 0.8932 -0.0238 175.21 0.0000
3 0.8523 0.0585 251.59 0.0000
4 0.8224 0.0526 323.45 0.0000
5 0.7907 0.0027 390.57 0.0000
6 0.7457 -0.2053 450.91 0.0000
7 0.6980 -0.1496 504.35 0.0000
8 0.6604 -0.0166 552.70 0.0000
9 0.6229 -0.1212 596.19 0.0000
10 0.5862 -0.0934 635.14 0.0000
...
40 -0.4325 -0.5224 1035.50 0.0000
Note
Slowly decaying AC values and a single large PAC spike at lag 1 are the hallmark of a non-stationary (unit root) process.
Logging variables is a frequently-recommended way to deal with the variance-non-stationarity of a variable.
In practice this doesn’t usually make much of a difference.
Of course if you have done this, you need to keep this in mind when interpreting your results.
| Model | ACF | PACF |
|---|---|---|
| White noise | All zero | All zero |
| MA(1): \theta > 0 | Single positive spike at lag 1 | Oscillating decay |
| MA(1): \theta < 0 | Single negative spike at lag 1 | Geometric decay |
| AR(1): \phi_1 > 0 | Direct geometric decay | Single spike at lag 1 |
| AR(1): \phi_1 < 0 | Oscillating decay | Single spike at lag 1 |
| AR(p) | Decays toward zero | Spikes through lag p; zero for s > p |
| ARMA(1,1) | Geometric decay after lag 1 | Oscillating decay after lag 1 |
| ARMA(p,q) | Decay after lag q | Decay after lag p |
Pure AR processes generally lead to declining ACFs and spikes in the PACF to indicate the ordering.
Pure MA processes show the opposite — a spike in the ACF and a declining pattern in the PACF.
Combinations (non-pure ARMA) are more difficult to diagnose.
A single spike in both the ACF and PACF can usually be modelled well by an MA at that lag.
Any declining pattern in the ACF indicates an AR process starting at the lag at which the initial spike appears.
Tip
Always keep in mind that we are moving from a position of ignorance. Significant results in the ACF and PACF may arise purely by chance.
A univariate time series can be characterised as ARIMA(p, d, q) where:
Once we know (p, d, q) for a particular series, we can write it out. For instance, if Y_t is ARIMA (1, 1, 1):
(1-L)Y_t = \phi(1-L)Y_{t-1} + \epsilon_t + \theta\epsilon_{t-1}
The Box-Jenkins approach requires that a variable be stationary and invertible.
This is why we start out determining whether or not the series is stationary.
A time series is said to be “invertible” if it can be represented by a finite-order or convergent autoregressive process.
Once we think we have the right (p, d, q) model, we test it with a filtering model.
One critique of R^2 is that it can only go up with the addition of new variables. Penalty-based criteria build in parsimony.
Adjusted R^2: \bar{R}^2 = 1 - \frac{n-1}{n-K}(1-R^2)
Akaike Information Criterion (AIC): \text{AIC}(K) = s_y^2(1-R^2)\,e^{2K/n}
Bayesian Information Criterion (BIC / Schwartz criterion): \text{BIC}(K) = s_y^2(1-R^2)\,n^{K/n}
For information criteria: smaller values indicate better-fitting, more parsimonious models. It is often suggested that AIC has better properties for forecasting models while BIC (generally) has a more extreme penalty for each parameter leading to more parsimonious explanations.
Starting from: Y_t - Y_{t-1} = \epsilon_t \;\;\Longrightarrow\;\; (1-L)Y_t = \epsilon_t \;\;\Longrightarrow\;\; (1-L)^1 Y_t = \epsilon_t
We can generalise by allowing d to take non-integer values:
(1-L)^d Y_t = \epsilon_t
| Value of d | Process type |
|---|---|
| -0.5 < d < 0 | Anti-persistent, stationary |
| 0 < d < 0.5 | Long-memoried, stationary (mean-reverting) |
| 0.5 \leq d < 1 | Long-memoried, non-stationary |
| d = 1 | Classic unit root (random walk) |
The ARFIMA model generalises the ARIMA model by allowing non-integer degrees of integration.
ARFIMA models provide a solution for the tendency to over-difference stationary series that exhibit long-run dependence.
In the standard ARIMA approach, d is constrained to integers (0, 1, 2, …). Many series exhibit too much dependence to be I(0) but are not I(1) — ARFIMA captures this middle ground.
Such series are said to be integrated of order d, denoted I(d).
The standard ARIMA model in lag-operator notation:
(1-L)Y_t = \frac{\theta(L)}{\phi(L)}\epsilon_t
The ARFIMA generalisation:
(1-L)^d Y_t = \frac{\theta(L)}{\phi(L)}\epsilon_t
where d can now range continuously from -0.5 to 1 (and beyond).
A simple special case — the Pure Fractional Noise model:
(1-L)^d Y_t = \epsilon_t, \quad \epsilon_t \sim N(0, \sigma^2)
Lebo, Walker, and Clarke (2000) provide evidence that many political time series are consistent with fractional integration.
Interventions or “impact assessments” are a popular extension of ARIMA modeling.
In the most basic models of interventions, we are testing a theory about the impact of a single event on our dependent variable of interest.
If we start out with a stationary dependent variable (Y_t), we can write an intervention model compactly as:
Y_t = f(I_{t-i}) + N_t
where:
Interventions in time series can be characterised by two dimensions:
Interventions can also be written generally as:
Y_t = \frac{\omega_0}{1 - \delta L}\, I_{t-i} + N_t
where \delta is known as the “adjustment parameter”.
For permanent intervention effects, I_t = 0 before the event and I_t = 1 thereafter.
E.g., I_t:
00000000011111111111
This formulation of I_t is also called a step function.
Abrupt permanent intervention (stationary Y_t, \delta = 0):
Y_t = \omega_0 I_{t-i} + N_t
This is a special case of the general model with \delta = 0.
If we keep the same step-function I_t but set 1 > \delta > 0, we have a gradual permanent intervention.
For temporary intervention effects, I_t = 0 before the event, I_t = 1 at the time of the event, and I_t = 0 thereafter.
E.g., I_t:
00000000010000000000
Abrupt temporary intervention (stationary Y_t, \delta = 0):
Y_t = \omega_0 I_{t-i} + N_t
This is again the general model with \delta = 0.
If we keep the same impulse I_t but set 0 < \delta < 1, we have a gradual temporary intervention.
Starting with the general formula for an ARIMA(0,0,0) Y_t (setting i = 0 and ignoring the noise term):
Y_t = \frac{\omega_0}{1 - \delta L}\, I_t
Multiply both sides by (1 - \delta L):
(1 - \delta L)\,Y_t = \omega_0 I_t
Expanding the lag operator:
Y_t - \delta Y_{t-1} = \omega_0 I_t
\boxed{Y_t = \delta\, Y_{t-1} + \omega_0 I_t}
This reveals Y_t as an AR(1) process driven by the intervention: each period’s value depends on the previous period scaled by \delta, plus the current intervention impulse.
Note: For \delta = 1 the step input induces a random walk (no long-run equilibrium).
It is also possible to model a gradual intervention while keeping \delta = 0:
Y_t = \omega_0 I_{t-i} + N_t
Instead of estimating \delta, we construct a ramp directly in I_t.
E.g., I_t:
0 0 0 .25 .5 .75 1 0 0 0
Caution: this approach imposes strong metric assumptions about the shape and timing of the transition — assumptions that are baked into the variable construction rather than estimated from the data.
Four canonical shapes (Enders, Applied Econometric Time Series, Figure 5.3):
| Panel | Label | Input | \delta | Description |
|---|---|---|---|---|
| (a) | Pure jump | Step | 0 | Abrupt permanent |
| (b) | Pulse | Impulse | 0 | Abrupt temporary |
| (c) | Gradually changing | Step | >0 | Gradual permanent |
| (d) | Prolonged pulse | Impulse | >0 | Gradual temporary |
The table maps directly onto the two-dimensional classification: dynamics (abrupt/gradual) × duration (temporary/permanent).
To model the impact of an intervention we first need a noise model for Y_t — but because we believe an intervention is present, this is potentially problematic.
Enders’ recommendation: Build the noise model using the longer of the two sub-series created by splitting the data at the point of the intervention I_{t-i}.
Alternative approach: Privilege the intervention variable —
There are multiple valid approaches; the key is not to let the intervention period contaminate the ARIMA identification step.
So far we have assumed a stationary Y_t. If Y_t is non-stationary and must be differenced, the intervention measure I_t must also be differenced.
Consider a permanent (step) intervention:
(1 - L)\,Y_t = \omega_0 I_t + \epsilon_t
Y_t - Y_{t-1} = \omega_0 I_t + \epsilon_t
Y_t = Y_{t-1} + \omega_0 I_t + \epsilon_t
The problem: we have differenced Y_t but not I_t. This introduces a random walk with drift when we actually want a random walk with a level shift. The intervention must be differenced to maintain consistency with the dependent variable transformation.
Box-Jenkins transfer function models extend the ARIMA framework to models with continuous independent variables.
In transfer function modelling, you either:
A central challenge: determining the appropriate lag at which to specify each independent variable.
Consider the model:
\Delta GOVPOPL_t = \omega_1 FALK_t + \omega_2 EMPXL_{t-i} + N_t
Different theories imply different values of i:
| i | Mechanism |
|---|---|
| 0 | “Market/me” — informational effect is immediate |
| 1 | Media disclosure of economic measures |
| 2 | Diffuse social network (“mates”) transmission |
Strategy 1 — Theory-driven
Use substantive knowledge about the underlying causal process to set i a priori.
Strategy 2 — Data-driven identification
Let the data guide lag specification through:
Pre-whitening controls for serial autocorrelation before examining cross-correlations:
Identify and fit the univariate ARIMA model for the independent variable X_t; save the residuals \hat{\alpha}_t.
Apply the same ARIMA filter to the dependent variable Y_t; save those residuals \hat{\beta}_t.
Cross-correlate \hat{\alpha}_t and \hat{\beta}_t to reveal the lag structure free of autocorrelation contamination.
Software note: RATS handles pre-whitening for AR (linear) specifications but not for MA (non-linear) processes. A second-best approach is to filter Y_t using its own univariate ARIMA and then cross-correlate with the pre-whitened X_t.
\Delta GOVPOPL_t = \omega_1 FALK_t + \omega_2 EMPXL_{t-i} + N_t
Suppose we expect the effect of unemployment (EMPXL) to be negative and believe i < 4 (effect materialises within one quarter).
Procedure:
This data-driven identification strategy is known as cross-correlation.
This complete analysis and the original article can be found on Box
The ARIMA approach is fundamentally inductive. The workflow involves the use of empirical values of ACFs and PACFs to engage in model selection. Dynamic models engage theory/structure to impose more stringent assumptions for producing estimates.
First, a result. Aitken Theorem
In a now-classic paper, Aitken generalized the Gauss-Markov theorem to the class of Generalized Least Squares estimators. It is important to note that these are GLS and not FGLS estimators. What is the difference? The two GLS estimators considered by Stimson are not strictly speaking GLS.
Definition \hat{\beta}_{GLS} = (\mathbf{X}^{\prime}\Omega^{-1}\mathbf{X})^{-1}\mathbf{X}^{\prime}\Omega^{-1}\mathbf{y} > Properties >
> (1) GLS is unbiased.
> (2) Consistent.
> (3) Asymptotically normal.
> (4) MV(L)UE
The variance/covariance matrix of the errors for a first-order autoregressive process is useful to derive.
The matrix is banded; observations separated by one point in time are correlated \rho. Period two is \rho^2; the corners are \rho^{T-1}. The diagonal is one.
What I have actually described is the correlation; the relevant autocovariances are actually defined by \frac{\sigma^{2}\rho^{s}}{1 - \rho^2} where s denotes the time period separation.
It is also straightforward to prove (tediously through induction) that this is invertible; it is square and the determinant is non-zero having assumed that |\rho < 1|.
\Phi = \sigma^{2}\Psi = \sigma^{2}_{e} \left(\begin{array}{ccccc}1 & \rho^{1} & \rho^{2} & \ldots & \rho^{T-1} \\ \rho^1 & 1 & \rho^1 & \ldots & \rho^{T-2} \\ \rho^{2} & \rho^1 & 1 & \ldots & \rho^{T-3} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ \rho^{T-1} & \rho^{T-2} & \rho^{T-3} & \ldots & 1 \end{array}\right)
given that e_{t} = \rho e_{t-1} + \nu_{t}. A Toeplitz form….
If the variance is stationary, we can rewrite, \sigma^{2}_{e} = \frac{\sigma^{2}_{\nu}}{1 - \rho^{2}}
A comment on characteristic roots….
We have the two key elements to implement this except that we do not know \rho; we will have to estimate it and estimates have uncertainty. But it is important to note this imposes exactly an AR(1). If the process is incorrectly specified, then the optimal properties do not follow. Indeed, the optimal properties also depend on an additional important feature.
We need to estimate things to replace unknown covariance structures and coverage will depend on properties of the estimators of these covariances.
Consistent estimators will work but there is euphemistically considerable variation in the class of consistent estimators.
Contrasting the Beck and Katz/White approach with the GLS approach is a valid difference in philosophies.1 One takes advantage of OLS and Basus Theorem; one goes full Aitken.
y_{t} = a_{1} y_{t-1} + \epsilon_{t}
is the simplest dynamic model but it cannot be estimated consistently, in general terms, in the presence of serial correlation. Why?
The key condition for unbiasedness is violated because \mathbb{E}(y_{t-1}\epsilon_{t}) \neq 0. OLS will not generally work.
A note on dynamic interpretation.
y_{t} = a_{1} y_{t-1} + \beta X_t + \epsilon_{t}
The problem is fitting and the key issue is white noise residuals post-estimation. But we have to assume a structure and implement it.
y_{t} = \alpha + \beta_{0} X_t + \beta_{1}x_{t-1} + \ldots + \epsilon_{t}
The impact of x occurs over multiple periods. It relies on theory, or perhaps analysis using information criteria/F [owing to quasi-nesting and missing data]. OLS is a fine solution to this problem but the search space of models is often large.
In response to this problem, we have structured distributed lag models; there are many such schemes.
Koyck/Geometric decay:
short run and long-run effects are parametrically identified y_t = \alpha + \beta(1-\lambda)\sum_{j=0}^{\infty}\lambda^{j}X_{t-j} + \epsilon
Almon (more arbitrary decay) y_{it} = \sum_{t_{A}=0}^{T_{F}} \rho_{t_{A}}x_{t - t_{A}} + \epsilon_{t} with coefficients that are ordinates of some general polynomial of degree T_{F} >> q. The \rho_{t_{A}} = \sum_{k=0}^{T_{F}} \gamma_{k}t^{k}.
y_{t} = \alpha + \gamma_{1}y_{t-1} + \beta_{0} X_t + \beta_{1}X_{t-1} + \beta_{2}X_{t-2} + \ldots + \epsilon_{t}
As recently as April of 2025, a paper appeared in the Journal of Politics advocating the use of ADL(2,2). The paper, by Kagalwala and Whitten, called The Answer was There All Along: Worry about the dynamics!. A previous argument was made for the ADL(1,1) a few years before. I have placed the 2025 paper in the Box.
Data analysis can quite yield models comparisons among competing dynamic structures. The key issue is that the analyst need divine the process; what is the relevant error process and what is the structure and timing of effects alongside the potential question of incremental adjustment. We need good theory for that.
Given such theory, we can take an equations as analysis approach, measure the variables, and derive reduced forms, and then recover parameter estimates deploying simultaneous equations methods. Very large such systems were a core part of early empirical macroeconomics. The failures of such systems led to the proposal of alternatives.
Chris Sims suggested a more flexible approach: the VAR.
The key insight is that this VAR is the reduced form to some more complicated as yet unspecified structural form.
But if the goal is to specify how variables related to one another and to use data to discover Granger causality and responses to impulse injected in the system.
library(forecast)
mdeaths
fdeaths
save(mdeaths, fdeaths, file = "./img/LungDeaths.RData")
Series: mdeaths, fdeaths
Model: VAR(3) w/ mean
Coefficients for mdeaths:
lag(mdeaths,1) lag(fdeaths,1) lag(mdeaths,2) lag(fdeaths,2)
0.6675 0.8074 0.3677 -1.4540
s.e. 0.3550 0.8347 0.3525 0.8088
lag(mdeaths,3) lag(fdeaths,3) constant
0.2606 -1.1214 538.7817
s.e. 0.3424 0.8143 137.1047
Coefficients for fdeaths:
lag(mdeaths,1) lag(fdeaths,1) lag(mdeaths,2) lag(fdeaths,2)
0.2138 0.4563 0.0937 -0.3984
s.e. 0.1460 0.3434 0.1450 0.3328
lag(mdeaths,3) lag(fdeaths,3) constant
0.0250 -0.315 202.0027
s.e. 0.1409 0.335 56.4065
Residual covariance matrix:
mdeaths fdeaths
mdeaths 58985.95 22747.94
fdeaths 22747.94 9983.95
log likelihood = -812.35
AIC = 1660.69 AICc = 1674.37 BIC = 1700.9
Series: mdeaths, fdeaths
Model: VAR(2) w/ mean
Coefficients for mdeaths:
lag(mdeaths,1) lag(fdeaths,1) lag(mdeaths,2) lag(fdeaths,2) constant
0.9610 0.3340 0.1149 -1.3379 443.8492
s.e. 0.3409 0.8252 0.3410 0.7922 124.4608
Coefficients for fdeaths:
lag(mdeaths,1) lag(fdeaths,1) lag(mdeaths,2) lag(fdeaths,2) constant
0.3391 0.2617 -0.0601 -0.2691 145.0546
s.e. 0.1450 0.3510 0.1450 0.3369 52.9324
Residual covariance matrix:
mdeaths fdeaths
mdeaths 62599.51 24942.79
fdeaths 24942.79 11322.70
log likelihood = -833.17
AIC = 1694.35 AICc = 1701.98 BIC = 1725.83
What happens if I shock one of the series; how does it work through the system?
The idea behind an impulse-response is core to counterfactual analysis with time series. What does our future world look like and what predictions arise from it and the model we have deployed?
Whether VARs or dynamic linear models or ADL models, these are key to interpreting a model in the real world.


ESSSSDA26-2J: One or More Time Series