Note that we don’t have exercise today. We will finish the part about counterfactual forecasting from last week’s exercise.
For Problem Set 2, the most important techniques are unit root test, model estimation, model selection, and counterfactual forecasting.
Although this Thursday Chris talked a little bit about panel data, I don’t want to go too far ahead. So I will move the panel data part to the next week.
I will return Problem Set 1 next Monday. Two thoughts here:
If we are not sure whether it is AR or MA process or how many order it should be, we can look at its longer-term behavior in ACF and PACF plots by specifying, for example, acf(ts, 100).
Besides decompose plot and PACF plot, we also need to overplot time series from each cycle to look at repetitive pattern.
Unit Root Test
Unit-root tests are the first step to study and analyze non-stationary time series.
Intuition: if time series is stationary, then regressing \(y_{t}-y_{t-1}\) on \(y_{t-1}\) should produce a negative coefficient. Why?
In a stationary series, knowing the past value of the series helps to predict the next period’s change. Positive shifts should be followed by negative shifts (mean reversion).
Augmented Dickey-Fuller test: null hypothesis of unit root.
Same with Phillips-Perron test, but differs in how the AR(\(p\)) time series is modeled w.r.t. lags, serial correlation, heteroskedasticity.
ADF Test
The ADF test regression can be expressed as follows:
\[\Delta y_t = \boldsymbol{\beta D_t} + \pi y_{t-1} + \sum^p_{j=1}\psi_j\Delta y_{t-j}+\epsilon_t\] where \(\boldsymbol{D}_t\) is a vector of deterministic trends. The \(p\) lagged difference terms, \(\Delta y_{t-j}\), are used to approximate the ARMA structure of the errors (consider the equivalence of AR and MA models), and the value of \(p\) is set so that the error is serially uncorrelated. The error term is assumed to be homoskedastic.
Under the null hypothesis (\(\rho - 1 = 0\) => unit root => non-stationary), \(\pi\) is set to 0 and is used to construct the ADF test statistic. If we don’t get a statistically significant result from ADF test, then the time series is likely to be stationary.
PP Test
The PP test regression can be expressed as follows:
The PP test ignores any serial correlation in the test regression. Instead, it corrects for serial correlation and heteroskedasticity in the errors \(u_t\) by directly modifying the test statistic, \(t_{\pi=0}\)
As we can see, both ADF test and PP test assume data has I(\(1\)) non-stationarity. But they also differ in how the time series is modeled. PP test allows for the residuals to be serially correlated.
Sometimes we may find that the results of ADF test and PP test are different. This is because some underlying temporal dependence is not incorporated into the assumed models, or because our sample is too small. For example, if there is a structural break (i.e. persistent arises due to large and infrequent shocks) in the time series, then ADF test tends to gives a non-stationary result.
If we suspect that our time series have particular forms of temporal dependence, we can use the following statistical tests: - higher-order differencing non-stationarity: Kwiatkowski-Phillips-Schmidt-Shin (KPSS) test - Structural break: Zivot-Andrews test, etc.
What is Non-stationary?
Recall that stationary time series have three properties:
Mean stationarity
mean does not depend on \(t\), constant over time
variance also does not depend on \(t\), constant over time
Covariance stationarity
covariance of \(y_t\) and \(y_{t-1}\) do not depend on \(t\)
does not depend on the actual time the covariance is computed
Ergodicity
sample moments coverge in probability to the population moments
sample mean and variance tend to be the same as entire series
Nonstationary procesess lack these properties. Note that we are abstracting from trends and other covariates.
What is Non-stationary?
Why do nonstationary processes matter?
ACF and PACF not defined since the covariances in the nominator depend on \(t\)
Spurious regression: we may detect strong correlations between nonstationary processes although they are really independent
Long run forecasts are very difficult since they do not tend toward any mean
What is Non-stationary?
Solutions?
Analyze nonstationary process using ARIMA (differencing)
effective at capturing short run changes
outcome is transformed to a difference
long-run predictions not feasible
Analyze nonstationary process using cointegration
effective at capturing long run relationships between nonstationary processes
outcome is left as a level
short-run and long-run predictions feasible
appropriate for analyzing multiple time series
Studying the Time Series
#Load packageslibrary(tseries) # For unit root testslibrary(forecast) # For decompose()library(lmtest) # For Breusch-Godfrey LM test of serial correlationlibrary(urca) # For estimating cointegration modelslibrary(simcf) # For counterfactual simulation via ldvsimev()library(MASS) # For mvrnorm()library(RColorBrewer) # For nice colorslibrary(Zelig) # For approval datalibrary(quantmod) # For creating lagssource("Lab5_data/TSplotHelper.R") # Helper function for counterfactual time series plots############################################################################Load data#US Presidential approval data (Bush, Monthly, 2/2001--6/2006)#Includes average oil price data ($/barrel?)data(approval)attach(approval)colnames(approval)
# plot presidential approvalplot(approve,type="l",ylab="Percent Approving",xlab="Time",main ="US Presidential Approval")lines(x=c(8,8),y=c(-1000,1000),col="red")lines(x=c(26,26),y=c(-1000,1000),col="blue")text("9/11",x =10, y =40, col="red",cex=0.7)text("Iraq \n War",x =28, y =40, col="blue",cex=0.7)
# plot average oil priceplot(avg.price,type="l",ylab="$ per Barrel",xlab="Time",main ="Average Price of Oil")lines(x=c(8,8),y=c(-1000,1000),col="red")lines(x=c(26,26),y=c(-1000,1000),col="blue")text("9/11",x =10, y =175, col="red",cex=0.7)text("Iraq \n War",x =28, y =175, col="blue",cex=0.7)
Unit Root Tests
Conflicting results from PP and ADF tests may suggest a structural break within the time series of presidential approval.
PP.test(approve)
Phillips-Perron Unit Root Test
data: approve
Dickey-Fuller = -2.8387, Truncation lag parameter = 3, p-value = 0.235
adf.test(approve)
Augmented Dickey-Fuller Test
data: approve
Dickey-Fuller = -3.9565, Lag order = 3, p-value = 0.01721
alternative hypothesis: stationary
PP.test(avg.price)
Phillips-Perron Unit Root Test
data: avg.price
Dickey-Fuller = -2.3318, Truncation lag parameter = 3, p-value = 0.4405
adf.test(avg.price)
Augmented Dickey-Fuller Test
data: avg.price
Dickey-Fuller = -3.0115, Lag order = 3, p-value = 0.1649
alternative hypothesis: stationary
Differencing the Time Series
Recall that we can use differencing to address the issue of nonstationarity.
Step 2: Generate Predicted Values and Prediction Intervals
approve_no911 <-predict(arima.res1a, newxreg = xhyp_no911) # predict after september 2001; it will return a list containing expected predictions and standard errors of predictions# construct a data.frame for plotpred_no911 <- tibble::tibble(mean =c(approve[1:7], approve_no911$pred),se =c(rep(0, 7), rep(approve_no911$se, 65-8+1)),lwr = mean -1.96*se, # we can use se here to construct lwr because we use tibble; you cannot do this if you use `data.frame()`upr = mean +1.96*se)DT::datatable(pred_no911)
# let's see how it goes compared to original time seriesplot(ts(approve, start =c(2001, 2), frequency =12))lines(ts(c(approve[7], approve_no911$pred), start =c(2001, 8), frequency =12), lty ="dashed")
Step 3: Simulate Expected Values and Confidence Intervals
# First we need the simulated parameterssims <-1e4simparams <-mvrnorm(n = sims, mu = arima.res1a$coef, Sigma = arima.res1a$var.coef)head(simparams)
Thus far, we have been examining a single time series with covariates. This assumes that there is no feedback between variables.
Yet, we may be interested in the relationship between two potentially nonstationary time series that influence each other.
The original idea behind cointegration is that two time series may be in equilibrium in the long run but in the short run the two series deviate from that equilibrium.
Cointegrated time series are two nonstationary time series that are causally connected and do not tend toward any particular level but tend toward each other.
Cointegration means that a specific combination of two nonstationary series may be stationary. We say that these two series are cointegrated and the vector that defines the stationary linear combination is called the cointegrating vector.
More on ECM Applications in political science:
Cointegration Analysis: Engle-Granger Two Step - Step 1
set.seed(98105)# Generate cointegrated datae1 <-rnorm(100)x <-cumsum(e1) # x is a cumulative sum of e1, so it must be non-stationarye2 <-rnorm(100)y <-0.6*x + e2 # define the relationship between x and y. They are both non-stationary#Run step 1 of the Engle-Granger two stepcoint.reg <-lm(y ~ x -1) #Estimate the cointegration vector by least squares with no constantcoint.err <-residuals(coint.reg) #This gives us the cotingeration vector.#Check for stationarity of the cointegration vectoradf.test(coint.err)$statistic |># get the test statistic of ADF testpunitroot(trend="nc") # punitroot() comes from urca package. it compute the cumulative probability of the ACF test statistic. "nc" means it is a test for regression without intercept (see ?urca::punitroot)
Dickey-Fuller
1.835502e-06
#Make the lag of the cointegration error termcoint.err.lag <- coint.err[1:(length(coint.err)-2)]#Make the difference of y and xdy <-diff(y)dx <-diff(x)#And their lagsdy.lag <- dy[1:(length(dy)-1)]dx.lag <- dx[1:(length(dx)-1)]#Delete the first dy, because we are missing lags for this obsdy <- dy[2:length(dy)]
Cointegration Analysis: Engle-Granger Two Step - Step 2
#Estimate an Error Correction Model with LSecm1 <-lm(dy ~ coint.err.lag + dy.lag + dx.lag)summary(ecm1)
Call:
lm(formula = dy ~ coint.err.lag + dy.lag + dx.lag)
Residuals:
Min 1Q Median 3Q Max
-3.05134 -0.64656 0.02462 0.76963 2.04272
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.08009 0.11676 0.686 0.494
coint.err.lag -0.89254 0.14294 -6.244 1.22e-08 ***
dy.lag -0.85110 0.11222 -7.584 2.36e-11 ***
dx.lag 0.56191 0.13088 4.293 4.28e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.152 on 94 degrees of freedom
Multiple R-squared: 0.4014, Adjusted R-squared: 0.3822
F-statistic: 21.01 on 3 and 94 DF, p-value: 1.688e-10
Cointegration Analysis: Johansen Estimator
#Alternatively, we can use the Johansen estimator#Create a matrix of the cointegrated variablescointvars <-cbind(y,x)# Perform cointegration testscoint.test1 <-ca.jo(cointvars,ecdet ="const",type="eigen",K=2,spec="longrun")summary(coint.test1)
######################
# Johansen-Procedure #
######################
Test type: maximal eigenvalue statistic (lambda max) , without linear trend and constant in cointegration
Eigenvalues (lambda):
[1] 3.520082e-01 2.399643e-02 1.460788e-17
Values of teststatistic and critical values of test:
test 10pct 5pct 1pct
r <= 1 | 2.38 7.52 9.24 12.97
r = 0 | 42.52 13.75 15.67 20.20
Eigenvectors, normalised to first column:
(These are the cointegration relations)
y.l2 x.l2 constant
y.l2 1.0000000 1.000000 1.0000000
x.l2 -0.5565254 -4.955866 0.4169868
constant -0.2979546 32.290892 32.4988636
Weights W:
(This is the loading matrix)
y.l2 x.l2 constant
y.d -0.92569359 0.003754100 -1.582064e-17
x.d 0.02021318 0.007824302 -4.297376e-18
ecm.res1 <-cajorls(coint.test1,r =1, # Cointegration rankreg.number =1) # which variable(s) to put on LHS#(column indexes of cointvars)summary(ecm.res1)
Length Class Mode
rlm 12 lm list
beta 3 -none- numeric
Cointegration Analysis
#For the presidential approval example, use an ECM equivalent to the#ARIMA(1,0,1) model that we chose earliercointvars <-cbind(approve,avg.price)ecm.test1 <-ca.jo(cointvars,ecdet ="const",type="eigen",K=2,spec="longrun",dumvar=cbind(sept.oct.2001,iraq.war))summary(ecm.test1)
######################
# Johansen-Procedure #
######################
Test type: maximal eigenvalue statistic (lambda max) , without linear trend and constant in cointegration
Eigenvalues (lambda):
[1] 2.391542e-01 1.367744e-01 1.387779e-16
Values of teststatistic and critical values of test:
test 10pct 5pct 1pct
r <= 1 | 9.27 7.52 9.24 12.97
r = 0 | 17.22 13.75 15.67 20.20
Eigenvectors, normalised to first column:
(These are the cointegration relations)
approve.l2 avg.price.l2 constant
approve.l2 1.0000000 1.0000000 1.00000
avg.price.l2 0.1535049 0.3382829 -1.05616
constant -76.0019182 -120.7778161 90.24200
Weights W:
(This is the loading matrix)
approve.l2 avg.price.l2 constant
approve.d -0.12619985 0.02231967 5.407866e-17
avg.price.d -0.02220281 -0.58771490 3.727850e-16
