Lecture 3: Time Series and Prediction
Have a first look at time series data
Discuss statistics measuring the persistence of univariate time series (autocovariance, autocorrelation)
Discuss various particular univariate processes
A key property: stationarity
Switch to bivariate time series and consider the danger of spurious regression
Consider more advanced time series models
Prediction in all aforementioned models
Material: Wooldridge Chapters 10, 11, 12, 18.5
Definition: Time Series Data
Time series data consists of observations of a variable or several variables over time.
Examples: Time Series Data
What Makes a Model “Dynamic”?
A common dynamic specification is the Autoregressive model:
\[ y_t = \alpha + \sum_{i=1}^{p} \beta_i y_{t-i} + \epsilon_t \]
Interpreting the coefficients (\(\beta_i\)) directly is often not the primary goal. Instead, we focus on key dynamic characteristics derived from them.
Key Questions to Answer When Interpreting Dynamic Models
A time trend is a variable, denoted as \(t\), that increases by one unit for each period (\(t = 1, 2, 3, \dots, T\)). It is included as an independent variable to capture sustained, long-term movements in the dependent variable.
A model with a linear time trend takes the form:
\[ y_t = \alpha + \beta t + \delta x_t + \epsilon_t \]
The coefficient \(\beta\) measures the average change in the dependent variable (\(y_t\)) per time period (e.g., per year or per quarter), holding all other included variables constant (ceteris paribus).
The trend term acts as a proxy for unobserved, slowly evolving factors like technological progress, population growth, or gradual shifts in consumer preferences.
Definition: Time Series (Mathematical)
A time series is a set of random variables indexed by time, denoted as \(\{Y_1, Y_2, \dots, Y_T\}\) or simply \(\{Y_t\}\).
Stochastic vs. Deterministic: We treat a time series as a stochastic process because we don’t know the future values with certainty. We can only talk about their probability distributions.
One Realization: The actual data we have (e.g., GDP from 1950-2024) is just one of many possible paths the process could have taken.
Our goal is to infer the properties of the underlying process from this single realization.
Just like any random variable, each point in a time series has statistical properties.
Crucially, in general, the mean and variance could be different at each point in time.
Definition: Autocovariance
The \(k\)-th Order Autocovariance (\(\gamma_k\)):
\[ \gamma(t, t-k) = Cov(Y_t, Y_{t-k}) = E[(Y_t - \mu_t)(Y_{t-k} - \mu_{t-k})] \]
The value of the autocovariance (\(\gamma_k\)) depends on the units of the variable \(Y\). This makes it hard to compare across different series.
Autocorrelation: standardize the autocovariance to get the autocorrelation, which is unit-free.
Definition: Autocorrelation
The \(k\)-th Order Autocorrelation (\(\rho_k\)):
\[ \rho_k = \frac{ Cov(Y_t, Y_{t-k}) }{\sqrt{Var(Y_t) Var(Y_{t-k})}} \]
If the series is stationary, this simplifies to:
\[ \rho_k = \frac{\gamma_k}{\gamma_0} \]
where \(\gamma_0\) is the variance, \(Var(Y_t)\).
Definition: White Noise Process
A white noise is a time series such that:
\[ Y_t = α + \rho Y_{t-1} + u_t \]
For analysis to be tractable, we often assume a time series is covariance stationary (or weakly stationary).
This means the three aforementioned statistical properties do not change over time. Three conditions must hold:
Definition: Stationarity
\[ y_t = \rho y_{t-1} + \epsilon_t \]
Directly testing \(\rho=1\) is problematic. Instead, Dickey and Fuller proposed a clever transformation. We subtract \(y_{t-1}\) from both sides of the AR(1) equation:
\[ y_t - y_{t-1} = \rho y_{t-1} - y_{t-1} + \epsilon_t \]
This simplifies to the classic Dickey-Fuller regression form:
\[ \Delta y_t = (\rho - 1) y_{t-1} + \epsilon_t \quad \implies \quad \Delta y_t = \gamma y_{t-1} + \epsilon_t \]
The test is now focused on the coefficient \(\gamma = (\rho - 1)\):
The simple DF test assumes the error term \(\epsilon_t\) is white noise. What if the true underlying process is more complex, like an AR(2), causing the errors in the simple DF regression to be serially correlated?
The Augmented Dickey-Fuller (ADF) Test adds as many lagged difference terms as necessary to eliminate serial correlation in the residuals of the test equation:
\[ \Delta y_t = \alpha + \beta t + \gamma y_{t-1} + \sum_{i=1}^{p} \delta_i \Delta y_{t-i} + \epsilon_t \]
# Import the function
from statsmodels.tsa.stattools import adfuller
# Assume 'y' is a pandas Series or NumPy array
# The most common test includes a constant ('c')
result = adfuller(y, regression='c', autolag='AIC')
# Print the key results
print(f'ADF Statistic: {result[0]}')
print(f'p-value: {result[1]}')
# result[1] is the p-value.
# To test with a constant and trend, use regression='ct'# Using the simpler 'tseries' package
# Install if needed: install.packages("tseries")
library(tseries)
# Assume 'y' is a time series object
# The default test includes a constant
result <- adf.test(y)
# Print the result
print(result)
# The p-value is directly reported.
# For a test with constant and trend using the more robust 'urca' package:
# library(urca)
# summary(ur.df(y, type = "trend", lags = ...))/* Assume the time series variable is named 'y' and you have tsset your data */
/* Test with a constant (the default) */
/* The 'lags(k)' option can be used to specify the number of lags */
dfuller y
/* Test with a constant and a trend */
dfuller y, trend
/* Stata reports the test statistic and critical values for 1%, 5%, and 10%.
You must compare the Test Statistic to the critical values.
To get a p-value, use the user-written command 'dfgls'. */The special case of non-stationarity in the AR(1) model occurs when \(\rho=1\):
\[ Y_t = Y_{t-1} + u_t \]
The best forecast for tomorrow’s value is today’s value.
Permanent Shocks: The impact of a shock \(u_t\) is permanent; it is carried forward in all future periods.
\[ Y_t = \alpha + \rho_1 Y_{t-1} + \rho_2 Y_{t-2} + \dots + \rho_p Y_{t-p} + u_t \]
In the MA model, the current value is a function of the current and past random shocks.
\[ Y_t = \mu + u_t + \theta u_{t-1} \]
\(\mu\): The mean of the series.
\(\theta\): The moving average coefficient.
The MA(1) process has a finite memory. A shock \(u_{t-1}\) affects \(Y_{t-1}\) and \(Y_t\), but has no direct effect on \(Y_{t+1}\).
# Generate sample data
set.seed(123)
data <- rnorm(100)
# Fit ARMA(1,1) model
arma_model <- arima(data, order = c(1, 0, 1)) # (p, d, q) where d=0 for ARMA
# Print summary
summary(arma_model)
## Length Class Mode
## coef 3 -none- numeric
## sigma2 1 -none- numeric
## var.coef 9 -none- numeric
## mask 3 -none- logical
## loglik 1 -none- numeric
## aic 1 -none- numeric
## arma 7 -none- numeric
## residuals 100 ts numeric
## call 3 -none- call
## series 1 -none- character
## code 1 -none- numeric
## n.cond 1 -none- numeric
## nobs 1 -none- numeric
## model 10 -none- listimport numpy as np
import statsmodels.api as sm
# Generate sample data
np.random.seed(123)
data = np.random.randn(100)
# Fit ARMA(1,1) model
arma_model = sm.tsa.ARIMA(data, order=(1,0,1))
arma_results = arma_model.fit()
# Print summary
print(arma_results.summary())
## SARIMAX Results
## ==============================================================================
## Dep. Variable: y No. Observations: 100
## Model: ARIMA(1, 0, 1) Log Likelihood -153.957
## Date: Wed, 29 Oct 2025 AIC 315.915
## Time: 12:56:11 BIC 326.335
## Sample: 0 HQIC 320.132
## - 100
## Covariance Type: opg
## ==============================================================================
## coef std err z P>|z| [0.025 0.975]
## ------------------------------------------------------------------------------
## const 0.0270 0.114 0.237 0.812 -0.196 0.250
## ar.L1 0.0035 16.295 0.000 1.000 -31.935 31.942
## ma.L1 0.0036 16.302 0.000 1.000 -31.949 31.956
## sigma2 1.2729 0.219 5.809 0.000 0.843 1.702
## ===================================================================================
## Ljung-Box (L1) (Q): 0.00 Jarque-Bera (JB): 1.74
## Prob(Q): 1.00 Prob(JB): 0.42
## Heteroskedasticity (H): 0.70 Skew: 0.03
## Prob(H) (two-sided): 0.31 Kurtosis: 2.36
## ===================================================================================
##
## Warnings:
## [1] Covariance matrix calculated using the outer product of gradients (complex-step).To build a dynamically complete model, we need a model whose residuals are white noise (i.e., free from autocorrelation), ensuring that we have captured all the systematic, predictable information in the data.
Step 1: Analyze the Variables (Pre-Modeling)
This step determines the correct form of the variables to include in our initial regression.
Definition: Breusch-Godfrey Test
The BG test assesses whether the error term, \(\epsilon_t\), follows an AR(p) process. You, the researcher, choose the number of lags (p) to test for. The assumed model for the error is: \[\epsilon_t = \rho_1 \epsilon_{t-1} + \rho_2 \epsilon_{t-2} + \dots + \rho_p \epsilon_{t-p} + v_t,\]
where \(v_t\) is a white noise term. \(H_0: \rho_1, \dots, \rho_p=0\), and \(H_A\) is that at least one \(\rho\) coefficient is non-zero: there is serial correlation up to some order \(p\).
Breusch-Godfrey Test: Procedure
The BG test is performed in three simple steps:
library(lmtest)
# Step 1: Simulate Data with Serially Correlated Errors
# For reproducibility of the random numbers
set.seed(42)
# Define simulation parameters
n <- 200 # Sample size
alpha <- 5 # True intercept
beta <- 2 # True slope for X
# *** Key step: Define the autocorrelation structure for the errors ***
# We will create an error term u_t = rho * u_{t-1} + e_t
# where e_t is white noise.
rho <- 0.8 # High positive autocorrelation coefficient
sigma_e <- 1.5 # Standard deviation of the white noise component
# Generate the independent variable X
X <- rnorm(n, mean = 20, sd = 10)
# Generate the serially correlated error term 'u'
u <- numeric(n)
e <- rnorm(n, mean = 0, sd = sigma_e) # White noise component
# The first error term has no preceding term
u[1] <- e[1]
# Generate the rest of the errors using the AR(1) formula
for (t in 2:n) {
u[t] <- rho * u[t-1] + e[t]
}
# Generate the dependent variable Y using the regression equation
# Y_t = alpha + beta*X_t + u_t
Y <- alpha + beta * X + u
# Combine into a data frame for easy use
sim_data <- data.frame(Y = Y, X = X)
# Step 2: Estimate the OLS Regression Model
# We estimate the model Y ~ X, pretending we don't know about the
# autocorrelation in the errors.
ols_model <- lm(Y ~ X, data = sim_data)
# Print the summary. The coefficients will likely be close to the true
# values, but the standard errors and p-values are unreliable.
print(summary(ols_model))
##
## Call:
## lm(formula = Y ~ X, data = sim_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.0295 -1.3534 -0.0239 1.3439 5.1120
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.52239 0.31093 17.76 <2e-16 ***
## X 1.97522 0.01414 139.70 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.944 on 198 degrees of freedom
## Multiple R-squared: 0.99, Adjusted R-squared: 0.9899
## F-statistic: 1.952e+04 on 1 and 198 DF, p-value: < 2.2e-16
#--------------------------------------------------------------------
# Step 3: Conduct the Breusch-Godfrey Test
#--------------------------------------------------------------------
# The bgtest() function from the 'lmtest' package performs the test.
# We need to specify the model and the order of autocorrelation to test for.
# Since we created an AR(1) process, we will test for order 1.
bg_test_order1 <- bgtest(ols_model, order = 1)
print(bg_test_order1)
##
## Breusch-Godfrey test for serial correlation of order up to 1
##
## data: ols_model
## LM test = 96.14, df = 1, p-value < 2.2e-16
# You can also test for higher orders. For example, testing up to order 3.
# This would test the joint null hypothesis that the coefficients on the first
# three lags of the residuals are all zero in the auxiliary regression.
bg_test_order3 <- bgtest(ols_model, order = 3)
print(bg_test_order3)
##
## Breusch-Godfrey test for serial correlation of order up to 3
##
## data: ols_model
## LM test = 96.794, df = 3, p-value < 2.2e-16import numpy as np
import pandas as pd
import statsmodels.api as sm
import statsmodels.formula.api as smf
from statsmodels.stats.diagnostic import acorr_breusch_godfrey
# Step 1: Simulate Data with Serially Correlated Errors
# For reproducibility of the random numbers
np.random.seed(42)
# Define simulation parameters
n = 200 # Sample size
alpha = 5 # True intercept
beta = 2 # True slope for X
# *** Key step: Define the autocorrelation structure for the errors ***
# We will create an error term u_t = rho * u_{t-1} + e_t
# where e_t is white noise.
rho = 0.8 # High positive autocorrelation coefficient
sigma_e = 1.5 # Standard deviation of the white noise component
# Generate the independent variable X
# In numpy, mean is 'loc' and standard deviation is 'scale'
X = np.random.normal(loc=20, scale=10, size=n)
# Generate the serially correlated error term 'u'
# Initialize arrays for the white noise and the final error term
e = np.random.normal(loc=0, scale=sigma_e, size=n)
u = np.zeros(n)
# The first error term has no preceding term (Python uses 0-based indexing)
u[0] = e[0]
# Generate the rest of the errors using the AR(1) formula
for t in range(1, n):
u[t] = rho * u[t-1] + e[t]
# Generate the dependent variable Y using the regression equation
# Y_t = alpha + beta*X_t + u_t
Y = alpha + beta * X + u
# Combine into a data frame for easy use
sim_data = pd.DataFrame({'Y': Y, 'X': X})
# Step 2: Estimate the OLS Regression Model
# We estimate the model Y ~ X, pretending we don't know about the
# autocorrelation in the errors.
# The formula syntax in statsmodels is identical to R's.
# We must call .fit() to perform the estimation.
ols_model = smf.ols('Y ~ X', data=sim_data).fit()
# Print the summary. The coefficients will likely be close to the true
# values, but the standard errors and p-values are unreliable.
print(ols_model.summary())
## OLS Regression Results
## ==============================================================================
## Dep. Variable: Y R-squared: 0.985
## Model: OLS Adj. R-squared: 0.985
## Method: Least Squares F-statistic: 1.273e+04
## Date: Wed, 29 Oct 2025 Prob (F-statistic): 1.24e-181
## Time: 12:56:11 Log-Likelihood: -452.92
## No. Observations: 200 AIC: 909.8
## Df Residuals: 198 BIC: 916.4
## Df Model: 1
## Covariance Type: nonrobust
## ==============================================================================
## coef std err t P>|t| [0.025 0.975]
## ------------------------------------------------------------------------------
## Intercept 5.3778 0.387 13.914 0.000 4.616 6.140
## X 2.0112 0.018 112.817 0.000 1.976 2.046
## ==============================================================================
## Omnibus: 4.086 Durbin-Watson: 0.452
## Prob(Omnibus): 0.130 Jarque-Bera (JB): 4.160
## Skew: 0.196 Prob(JB): 0.125
## Kurtosis: 3.587 Cond. No. 50.7
## ==============================================================================
##
## Notes:
## [1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
# Step 3: Conduct the Breusch-Godfrey Test
# The acorr_breusch_godfrey() function performs the test.
# It takes the fitted model results and the order of autocorrelation (nlags).
# It returns (lm_statistic, lm_pvalue, f_statistic, f_pvalue)
# Test for order 1, as we created an AR(1) process.
bg_test_order1 = acorr_breusch_godfrey(ols_model, nlags=1)
print("\n" + "="*50)
##
## ==================================================
print("Breusch-Godfrey Test (Order 1)")
## Breusch-Godfrey Test (Order 1)
print(f" LM Statistic: {bg_test_order1[0]:.4f}")
## LM Statistic: 119.7827
print(f" p-value (LM): {bg_test_order1[1]:.4f}")
## p-value (LM): 0.0000
# A very low p-value indicates rejection of the null hypothesis (no serial correlation).
# You can also test for higher orders. For example, testing up to order 3.
# This tests the joint null hypothesis that the coefficients on the first
# three lags of the residuals are all zero in the auxiliary regression.
bg_test_order3 = acorr_breusch_godfrey(ols_model, nlags=3)
print("\n" + "="*50)
##
## ==================================================
print("Breusch-Godfrey Test (Order 3)")
## Breusch-Godfrey Test (Order 3)
print(f" LM Statistic: {bg_test_order3[0]:.4f}")
## LM Statistic: 119.9787
print(f" p-value (LM): {bg_test_order3[1]:.4f}")
## p-value (LM): 0.0000
print("="*50)
## ==================================================* Housekeeping: Clear memory and stop the output from pausing
clear all
set more off
* Step 1: Simulate Data with Serially Correlated Errors
set seed 42
* Define simulation parameters using local macros
local n = 200 // Sample size
local alpha = 5 // True intercept
local beta = 2 // True slope for X
local rho = 0.8 // High positive autocorrelation coefficient
local sigma_e = 1.5 // Standard deviation of the white noise component
* Set the number of observations for the dataset
set obs `n'
* Generate an index variable 't' to declare this as time-series data
gen t = _n
tsset t
* Generate the independent variable X
gen X = rnormal(20, 10)
* Generate the serially correlated error term 'u'
* First, create the white noise component 'e'
gen e = rnormal(0, `sigma_e')
* Now, generate 'u' using the recursive AR(1) formula
* u_t = rho * u_{t-1} + e_t
* We initialize u as missing, set the first value, then generate the rest
gen u = .
replace u = e in 1
replace u = `rho' * L.u + e in 2/l // L. is the lag operator, 2/l means "from obs 2 to the last"
* Generate the dependent variable Y using the regression equation
* Y_t = alpha + beta*X_t + u_t
gen Y = `alpha' + `beta'*X + u
* Step 2: Estimate the OLS Regression Model
* We estimate the model Y ~ X, pretending we don't know about the
* autocorrelation in the errors. The output is printed automatically.
regress Y X
* Step 3: Conduct the Breusch-Godfrey Test
* The `estat bgodfrey` command is a post-estimation command that must be
* run immediately after the regression.
* Test for order 1, as we created an AR(1) process.
estat bgodfrey, lags(1)
* You can also test for higher orders. For example, testing up to order 3.
* This tests the joint null hypothesis that the coefficients on the first
* three lags of the residuals are all zero in the auxiliary regression.
estat bgodfrey, lags(3)The starting point for a relationship between time series is often the linear model.
The model form is identical to its cross-sectional counterpart:
\[ Y_t = \alpha + \beta X_t + \epsilon_t \]
This model posits a contemporaneous relationship between X and Y.
Question: Can we estimate this with OLS and trust the results?
For OLS to be a good estimator, the standard assumptions must hold. The most problematic one for time series data is the assumption of finite variance:
\[ \text{Var}(X_t) = \sigma^2 \]
Let’s look at the most common non-stationary process: the random walk.
\[ X_t = X_{t-1} + u_t \]
Theorem: Variance of a Random Walk
Suppose \(X_0=0\). Let us find the variance of \(X_t\).
We know that \(X_1 = u_i\). \(X_2 = X_1 + u_2 = u_1 + u_2\), etc.
Hence, \(X_t = \sum u_i\), and \(Var(X_t) = Var(u_1 + u_2 + \dots + u_t)\).
Since the \(u_i\) are uncorrelated (white noise), the variance of the sum is the sum of the variances: \[ \begin{align} Var(X_t) &= \sigma^2_u + \dots + \sigma^2_u \text{ ( t times)} \\ &= t \times \sigma^2_u \end{align} \]
The variance depends on time \(t\) and grows without bound. This is a clear violation of stationarity.
If the static model \(Y_t = \alpha + \beta X_t + \epsilon_t\) is so dangerous, what is the alternative?
Differencing the Data: If two variables are non-stationary, their first differences (\(\Delta Y_t\) and \(\Delta X_t\)) may be stationary.
Dynamic models: We need models that explicitly account for dynamics—the idea that a change in X can affect Y both today and in the future.
Definition: Distributed Lag Model
A finite distributed lag model of order \(q\) is:
\[Y_t = \alpha + \beta_0 X_t + \beta_1 X_{t-1} + \dots + \beta_q X_{t-q} + u_t\] where:
\(\beta_0\): The impact multiplier - the immediate effect of a one-unit change in \(X_t\) on \(Y_t\).
\(\beta_s (s > 0)\): The dynamic multipliers - the effect of a one-unit change in \(X_{t-s}\) on \(Y_t\).
\[Y_t = \alpha + \beta_0 X_t + \beta_1 X_{t-1} + \dots + \beta_q X_{t-q} + u_t\]
Definition: ARDL Model
An ARDL(p,q) model is defined as follows:
\[Y_t = \alpha + \sum(ρ_j \cdot Y_{t-j}) \text{ [from j=1 to p]} + \sum(\beta_j \cdot X_{t-m}) \text{ [from m=0 to q]} + u_t\]
The individual coefficients represent direct, short-term impacts before feedback effects occur.
Impact Multiplier (\(\beta_0\)): The immediate effect of a change in \(X_t\) on \(Y_t\).
Delay Multipliers (\(\beta_j\)): The effect of past changes in \(X\) on current \(Y\).
Autoregressive Coefficients (\(\phi_i\)): Measure the persistence or inertia in the process.
The LRM shows the total effect on Y after a permanent change in X has fully worked through the system.
Derivation (via Equilibrium Condition):
The Long-Run Multiplier (\(\theta\)):
Interpretation: A permanent one-unit increase in \(X\) is associated with a total long-run change of \(\theta\) units in \(Y\), after all dynamic adjustments are complete.
The ARDL model powerfully separates the immediate impacts (\(\beta_j\)) from the total equilibrium effect (\(\theta\)), which accounts for the feedback from Y’s own past values.
ARDL package is the standard for estimating these models.ardl command (in newer versions) or a widely used user-written one (ssc install ardl). It automates lag selection and estimation.statsmodels library has incorporated a dedicated ARDL class, which makes estimation much more straightforward than the previous manual method of creating lagged variables.library(ARDL)
library(urca) # Contains the example 'denmark' dataset
# Load the data
data(denmark)
# The variables are LRM, LRY, IBO, IDE (inflation)
# Find the optimal lag order automatically
# We model LRM as a function of LRY and IBO, with a max lag of 4 for each
auto_model <- auto_ardl(LRM ~ LRY + IBO, data = denmark, max_order = 4)
# Print the selected orders
print(auto_model$best_order)
## LRM LRY IBO
## 3 1 0
# Estimate the selected ARDL model
# Let's say auto_ardl selected ARDL(3, 1, 4)
ardl_model <- ardl(LRM ~ LRY + IBO, data = denmark, order = c(3, 1, 4))
# Get a detailed summary of the model
# This output neatly separates short-run and long-run coefficients
print(summary(ardl_model))
##
## Time series regression with "ts" data:
## Start = 5, End = 55
##
## Call:
## dynlm::dynlm(formula = full_formula, data = data, start = start,
## end = end)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.033673 -0.012617 -0.002257 0.008639 0.073178
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.83728 0.68121 2.697 0.010189 *
## L(LRM, 1) 0.46368 0.14626 3.170 0.002920 **
## L(LRM, 2) 0.56593 0.14483 3.908 0.000351 ***
## L(LRM, 3) -0.28115 0.11100 -2.533 0.015338 *
## LRY 0.57978 0.14671 3.952 0.000307 ***
## L(LRY, 1) -0.36536 0.15667 -2.332 0.024817 *
## IBO -1.14595 0.35807 -3.200 0.002688 **
## L(IBO, 1) 0.01182 0.61215 0.019 0.984693
## L(IBO, 2) 0.59817 0.63372 0.944 0.350890
## L(IBO, 3) -0.81318 0.56768 -1.432 0.159783
## L(IBO, 4) 0.37907 0.37942 0.999 0.323763
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.02095 on 40 degrees of freedom
## Multiple R-squared: 0.9848, Adjusted R-squared: 0.981
## F-statistic: 259.3 on 10 and 40 DF, p-value: < 2.2e-16import pandas as pd
import statsmodels.api as sm
from statsmodels.tsa.api import ARDL, ardl_select_order
# Load the same 'lutkepohl2' dataset from statsmodels datasets
data = r.denmark
data.index = pd.to_datetime(data['ENTRY'], format='%Y:%m').dt.to_period('Q')
# Select the variables of interest
data_subset = data[['LRM', 'LRY', 'IBO']]
# Assign endogenous (Y) and exogenous (X) variables
endog = data_subset['LRM']
exog = data_subset[['LRY', 'IBO']]
# Find the optimal lag order
# maxlag is for the AR (dependent) part, maxorder is for the DL (independent) part
# aic is the information criterion to use for selection
sel_order = ardl_select_order(
endog, maxlag=4, exog=exog, maxorder=4, ic='aic', trend='c'
)
# Instantiate the ARDL model with the selected order
# sel_order.ardl_order gives the optimal (p, q1, q2, ...)
ardl_result = sel_order.model.fit()
# Print the summary of the results
# This shows the short-run coefficients
print(ardl_result.summary())
## ARDL Model Results
## ==============================================================================
## Dep. Variable: LRM No. Observations: 55
## Model: ARDL(3, 1, 0) Log Likelihood 132.589
## Method: Conditional MLE S.D. of innovations 0.019
## Date: Wed, 29 Oct 2025 AIC -249.178
## Time: 12:56:12 BIC -233.568
## Sample: 06-30-1974 HQIC -243.194
## - 03-31-1987
## ==============================================================================
## coef std err z P>|z| [0.025 0.975]
## ------------------------------------------------------------------------------
## const 1.9372 0.364 5.326 0.000 1.205 2.670
## LRM.L1 0.4330 0.118 3.661 0.001 0.195 0.671
## LRM.L2 0.5434 0.129 4.213 0.000 0.284 0.803
## LRM.L3 -0.2511 0.093 -2.696 0.010 -0.439 -0.063
## LRY.L0 0.6018 0.132 4.574 0.000 0.337 0.867
## LRY.L1 -0.3562 0.141 -2.521 0.015 -0.641 -0.072
## IBO.L0 -1.0497 0.178 -5.883 0.000 -1.409 -0.690
## ==============================================================================
params = ardl_result.params
denominator = 1 - params.filter(like='LRM.L').sum()
# 3. Calculate the sum of coefficients for each independent variable and divide
long_run_effects = pd.Series({
'LRY': params.filter(like='LRY.L').sum() / denominator,
'IBO': params.filter(like='IBO.L').sum() / denominator,
'const': params['const'] / denominator # Optional: long-run intercept
})
print(long_run_effects)
## LRY 0.893800
## IBO -3.821235
## const 7.052197
## dtype: float64* Load the data (same as the 'denmark' dataset in R)
webuse lutkepohl2, clear
* Declare the data as time series
tsset qtr
* Estimate the ARDL model, letting Stata select the best lags up to a max of 4.
* The dependent variable is 'ln_rm' (log real money).
* The independent variables are 'ln_inc' (log income) and 'ln_consump' (log consumption).
* We'll use the variables from the Stata dataset.
ardl ln_rm ln_inc, maxlags(4)
* After estimation, Stata shows the optimal lags and the short-run results.
* To see the long-run coefficients, use the post-estimation command `estat ec`.
estat ecDefinition: Forecasting
To predict future values of a time series based on its own past and potentially other related variables.
We use the information available up to time \(t\) to form the best possible guess about the value of the series at a future time \(t+h\).
We want to find the conditional expectation of the future value, given our current information.
\[ \text{Forecast} = E[Y_{t+h} | I_t] \]
Definition: Adjusted \(R^2\)
\[ R^2_{adj} = 1 - \frac{RSS / (n - k - 1)}{TSS / (n - 1)} \]
Where \(n\) = number of observations, and \(k\) = number of predictors (independent variables)
This can also be expressed in terms of the standard \(R^2\):
\[ R^2_{adj} = 1 - (1 - R^2) \frac{n - 1}{n - k - 1} \]
Definition: AIC
\[AIC = 2k - 2ln(L)\]
Where:
Note: For OLS models assuming normal errors, the AIC can be expressed using the RSS: \(AIC = n \cdot \ln(RSS/n) + 2k\)
The AIC provides a systematic way to select the optimal number of lags. The general procedure is as follows:
# Load a sample time series dataset
data(sunspot.year)
# Fit an AR model and let the function select the best lag order based on AIC
ar_model <- ar(sunspot.year)
# The optimal number of lags is stored in the 'order' component of the model object
optimal_lags <- ar_model$order
# Print the optimal number of lags
print(paste("Optimal number of lags based on AIC:", optimal_lags))
## [1] "Optimal number of lags based on AIC: 9"
# You can also view the AIC values for different lag orders
print(ar_model$aic)
## 0 1 2 3 4 5 6
## 500.451015 188.271695 37.689005 31.834674 33.427738 35.353944 28.912532
## 7 8 9 10 11 12 13
## 23.654977 9.099444 0.000000 1.973243 3.377547 5.376389 7.146043
## 14 15 16 17 18 19 20
## 8.037731 7.970507 9.526070 5.107497 6.785910 8.666711 10.661378
## 21 22 23 24
## 10.633330 12.467419 13.087979 14.552615import pandas as pd
from statsmodels.tsa.ar_model import ar_select_order
import statsmodels.api as sm
# Load a sample time series dataset (e.g., Sunspots data)
# This data is available in statsmodels
data = sm.datasets.sunspots.load_pandas().data
series = data['SUNACTIVITY']
# Use ar_select_order to find the optimal lag
# We specify a maximum lag to consider and the information criterion to use ('aic')
selection = ar_select_order(series, maxlag=20, ic='aic')
# The selected lags are stored in the 'ar_lags' attribute
optimal_lags = selection.ar_lags
# Print the optimal number of lags
print(f"Optimal number of lags based on AIC: {optimal_lags}")
## Optimal number of lags based on AIC: [1, 2, 3, 4, 5, 6, 7, 8, 9]* Load a sample time series dataset (Sunspot data)
* webuse is Stata's command to load data from the internet
webuse sunspot
* Declare the data as a time series
* The 'year' variable indicates the time period
tsset year
* Use varsoc to get the selection criteria for lags 1 through 20
* Stata will automatically mark the optimal lag for each criterion with a '*'
varsoc sunspot, maxlag(20)Definition: BIC
\[ BIC = -2 \ln(\hat{L}) + k' \ln(n) \]
Where:
Note: For OLS models assuming normal errors, the BIC can be expressed using the RSS: \(BIC = n \cdot \ln(RSS/n) + k' \ln(n)\)
Definition: Forecast Error
The forecast error is the difference between an actual observation and a forecast.
\[ e_t = Y_t - \hat{Y}_t \]
Sometimes the scale of the data matters. An error of 10 is huge for a series with an average of 5, but tiny for a series with an average of 5,000.
Scale-Free Metric:
In-Sample vs. Out-of-Sample Evaluation:
Let’s start predicting with a simple AR(p) model.
\[ Y_t = c + \phi_1 Y_{t-1} + \phi_2 Y_{t-2} + \dots + \phi_p Y_{t-p} + \varepsilon_t \]
Produce forecasts by iterative substitution:
For stationary AR processes, the long-run forecast (as \(h \to \infty\)) will converge to the unconditional mean of the process: \(E[Y] = \frac{c}{1 - \sum \phi_i}\).
Example: GDP Growth
Let’s assume we are modeling a country’s quarterly GDP growth rate, \(g_t\). After analyzing the data, we have estimated the following AR(1) model:
\[ g_t = 0.20 + 0.70 g_{t-1} + \epsilon_t \]
Where:
We are at the end of Quarter 4, 2024 (let’s call this period T). The last observed GDP growth rate was 1.5%.
Our goal is to forecast GDP growth for Q1 2025 (period T+1) and Q2 2025 (period T+2).
Example: GDP Growth (Continued)
A. Forecast for Q1 2025 (1-Step-Ahead Forecast)
To forecast for period T+1, we use the actual, known data from period T. The forecast, denoted with a “hat”, is:
\[ \hat{g}_{T+1} = c + \phi_1 g_T \]
Plugging in our known values:
\[ \hat{g}_{T+1} = 0.20 + 0.70 \times (1.5) \]
\[ \hat{g}_{T+1} = 0.20 + 1.05 = 1.25 \]
Our forecast for Q1 2025 growth is 1.25%.
Example: GDP Growth (Continued)
B. Forecast for Q2 2025 (2-Step-Ahead Forecast)
To forecast for period T+2, we need the value from the previous period (T+1). Since we are at time T, we do not know the actual value of \(g_{T+1}\). Therefore, we must use our own forecast, \(\hat{g}_{T+1}\), as the input.
\[ \hat{g}_{T+2} = c + \phi_1 \hat{g}_{T+1} \]
Plugging in the forecast we just calculated:
\[ \hat{g}_{T+2} = 0.20 + 0.70 \times (1.25) \]
\[ \hat{g}_{T+2} = 0.20 + 0.875 = 1.075 \]
Our forecast for Q2 2025 growth is 1.075%.
Example: GDP Growth (Continued)
| Period | Forecast (\(\hat{g}\)) | Actual (\(g\)) |
|---|---|---|
| Q1 2025 | 1.25% | 1.40% |
| Q2 2025 | 1.075% | 1.10% |
Step 1: Calculate the Forecast Errors (\(e_t = g_t - \hat{g}_t\))
Example: GDP Growth (Continued)
Step 2: Square the Errors and Calculate the Mean (MSE)
The Mean Squared Error (MSE) is the average of the squared errors.
\[ MSE = \frac{1}{n} \sum_{i=1}^{n} e_i^2 \]
\[ MSE = \frac{1}{2} \left( (0.15)^2 + (0.025)^2 \right) \]
\[ MSE = \frac{1}{2} \left( 0.0225 + 0.000625 \right) = \frac{0.023125}{2} = 0.0115625 \]
Example: GDP Growth (Continued)
Step 3: Take the Square Root of the MSE (RMSE)
The Root Mean Squared Error (RMSE) brings the metric back into the original units of the data (in this case, percentage points).
\[ RMSE = \sqrt{MSE} \]
\[ RMSE = \sqrt{0.0115625} \approx 0.1075 \]
Example: GDP Growth (Continued)
The RMSE for our two-period forecast is 0.1075%.
This value can be interpreted as the typical, or average, magnitude of our forecast error.
It means that, over this forecast horizon, our AR(1) model’s predictions were, on average, off by about 0.11 percentage points.
This single number provides a concise summary of the overall predictive accuracy of the model for this specific forecast period.
Now for an MA(q) model. The “memory” of these models is very different.
\[ Y_t = \mu + \varepsilon_t + \theta_1 \varepsilon_{t-1} + \dots + \theta_q \varepsilon_{t-q} \]
For any forecast horizon \(h\) greater than the order \(q\) (\(h > q\)), all the error terms in the equation will be in the future.
An ARMA model combines the features of both AR and MA processes.
\[ Y_t = c + \phi_1 Y_{t-1} + \dots + \phi_p Y_{t-p} + \varepsilon_t + \theta_1 \varepsilon_{t-1} + \dots + \theta_q \varepsilon_{t-q} \]
Forecasting combines the two previous methods:
Example: ARMA(1,1) Forecast
The AR component gives the forecast a persistent, decaying memory, while the MA component captures the effect of short-term, finite shocks.
Why limit ourselves to just the past of Y? Other variables might help predict it.
\[ Y_t = c + \sum_{i=1}^{p} \phi_i Y_{t-i} + \sum_{j=0}^{q} \beta_j X_{t-j} + \varepsilon_t \]
The process is similar to a standard AR forecast, but with an added component.
To forecast \(Y_{t+1}\), the model requires the value of \(X_{t+1}\) if \(\beta_0 \neq 0\)!
To forecast \(Y_{t+h}\), you need values for \(X_{t+h}, X_{t+h-1}, \dots, X_{t+1}\).
Solutions:
The Foundation is Conditional Expectation: all standard forecasting methods are an attempt to calculate \(E[Y_{t+h} | I_t]\).
Test your forecasts on out-of-sample data using metrics like MAE, RMSE, or MAPE to understand their real-world performance.
Model Structure dictates forecast behavior:
This lecture has been an introduction. More advanced topics include:
Empirical Economics: Lecture 3 - Time Series and Prediction