\documentclass[notitlepage]{article}
\usepackage{../../Math545}
\usepackage{listings}
\usepackage{numprint}
\usepackage{enumitem}
\usepackage{bbm}
\usepackage{bm}
\usepackage{verbatim}

\usepackage{amsfonts,amsmath,amssymb}

\newtheorem{alg}{Algorithm}
\newtheorem{ex}{Example}
\newtheorem{note}{Note}
\newenvironment{proof}{\par\noindent{\normalfont{\bf{Proof
}}}}
{\hfill \small$\blacksquare$\\[2mm]}


\setlength{\parindent}{0in}

\def\E{\mathbb{E}}
\def\Var{\text{Var}}
\def\Cov{\text{Cov}}
\def\Corr{\text{Corr}}

\def\bW{\mathbf{W}}

\newcommand{\bs}[1]{\bm{#1}}
\newcommand{\Rho}{\mathrm{P}}
\newcommand{\cif}{\text{ if }}
\newcommand{\Ra}{\Longrightarrow}
\newcommand{\ra}{\longrightarrow}

\def\bphi{\boldsymbol \phi}


\lstloadlanguages{R}

\definecolor{keywordcolor}{rgb}{0,0.6,0.6}
\definecolor{delimcolor}{rgb}{0.461,0.039,0.102}
\definecolor{Rcommentcolor}{rgb}{0.101,0.043,0.432}

\lstdefinestyle{Rsettings}{
  basicstyle=\ttfamily,
  breaklines=true,
  showstringspaces=false,
  keywords={if, else, function, theFunction, tmp}, % Write as many keywords
  otherkeywords={},
  commentstyle=\itshape\color{Rcommentcolor},
  keywordstyle=\color{keywordcolor},
  moredelim=[s][\color{delimcolor}]{"}{"},
}

\lstset{basicstyle=\ttfamily, numbers=none, literate={~} {$\sim$}{2}}

\def\ctitle{MATH 545: Introduction to Time Series Analysis}
\def\ttitle{Winter 2025}
\def\stitle{Cointegration, Error Correction and Granger Causality}

\begin{document}

<<setup, cache=FALSE, echo=FALSE>>=
library(knitr)
# global chunk options
opts_chunk$set(cache=TRUE, autodep=TRUE)
options(scipen=999)
options(repos=c(CRAN="https://cloud.r-project.org/"))
@

\maketitle

\textbf{Example:} Consider the VARMA$(1,1)$ model
\[
\bPhi(B) \bX_t =  \bTh (B) \bZ_t
\]
with
\begin{eqnarray*}
\bPhi(B) = \bPhi_1(B) = \left[ \begin{array}{cc} 1- 0.5 B & B \\ 0.25 B & 1-0.5 B \end{array} \right] \qquad \qquad
\bTh(B) = \bTh_1(B) = \left[ \begin{array}{cc} 1- 0.2 B & 0.4 B \\ 0.1 B & 1-0.2B \end{array} \right]
\end{eqnarray*}
with $\Var[\bZ_t] = \Sigma$.  That is
\[
\left[ \begin{array}{c} X_{1,t} \\ X_{2,t} \end{array} \right] = 
\left[ \begin{array}{rr} 0.5 & -1 \\ -0.25  & 0.5 \end{array} \right]
\left[ \begin{array}{c} X_{1,t-1} \\ X_{2,t-1} \end{array} \right] + \left[ \begin{array}{c} Z_{1,t} \\ Z_{2,t} \end{array} \right] + \left[ \begin{array}{rr} -0.2 & 0.4 \\ 0.1  & -0.2 \end{array} \right] \left[ \begin{array}{c} Z_{1,t-1} \\ Z_{2,t-1} \end{array} \right]
\]
or
\[
\bX_t = \bPhi \bX_{t-1} + \bZ_t + \bTh \bZ_{t-1}.
\]
say.  Then
\[
\left[ \begin{array}{cc} 1- 0.5 B & B \\ 0.25 B & 1-0.5 B \end{array} \right] \bX_t =
\left[ \begin{array}{cc} 1- 0.2 B & 0.4 B \\ 0.1 B & 1-0.2B \end{array} \right] \bZ_t
\]
This process is nonstationary, as we can discover by finding the eigenvalues of the $\bPhi$ matrix:
<<Ex1>>=
Phi<-matrix(c(0.5,-1,-0.25,0.5),2,2,byrow=T)
eigen(Phi)$values
@
Filtering on both sides by
\[
\left[ \begin{array}{cc} 1- 0.5 B & -B \\ -0.25 B & 1-0.5 B \end{array} \right]
\]
we obtain on the left and right hand sides respectively
\[
\left[ \begin{array}{cc} 1- 0.5 B & -B \\ -0.25 B & 1-0.5 B \end{array} \right] \left[ \begin{array}{cc} 1- 0.5 B & B\\ 0.25 B & 1-0.5 B \end{array} \right] = \left[ \begin{array}{cc} 1- B &  1\\ 0 & 1-B \end{array} \right]
\]
and
\[
\left[ \begin{array}{cc} 1- 0.5 B & -B \\ -0.25 B & 1-0.5 B \end{array} \right] \left[ \begin{array}{cc} 1- 0.2 B & 0.4 B\\ 0.1 B & 1-0.2 B \end{array} \right] = \left[ \begin{array}{cc} 1- 0.7 B & -0.6B \\ -0.15 B & 1-0.7 B \end{array} \right].
\]
Therefore
\[
\left[ \begin{array}{cc} 1- B & 0 \\ 0 & 1-B \end{array} \right] \bX_t
= \left[ \begin{array}{cc} 1- 0.7 B & -0.6B \\ -0.15 B & 1-0.7 B \end{array} \right] \bZ_t
\]
That is
\[
\left[ \begin{array}{l} (1- B)X_{1,t} \\  (1-B)X_{2,t}\end{array} \right]
= \left[ \begin{array}{l} Z_{1,t} - 0.7 Z_{1,t-1} - 0.6 Z_{2,t-1} \\ Z_{2,t} -0.15 Z_{1,t-1} - 0.7 Z_{2,t-1} \end{array} \right]
\]
Procesess $X_{1,t}$ and $X_{2,t}$ have unit roots, specifically
\[
X_{1,t},X_{2,t} \sim \text{ARIMA}(0,1,1).
\]
For example, with
\[
\Sigma = \left[ \begin{array}{cc} 0.50 & 0.0 \\ 0.0 & 0.50 \end{array} \right]
\]
<<Ex2,fig.height=4.5,fig.width=7>>=
set.seed(132)
library(mvnfast)
n<-1000
Sigma<-matrix(c(0.5,0,0,0.5),nrow=2,byrow=T)
Phi<-matrix(c(0.5,-1,-0.25,0.5),2,2,byrow=T)
Theta<-matrix(c(-0.2,0.4,0.1,-0.2),2,2,byrow=T)
Zt<-rmvn(n,c(0,0),Sigma)
Xt<-matrix(0,nrow=n,ncol=2)
for(i in 2:n){
    Xt[i,]<-Phi %*% Xt[i-1,] + Zt[i,] + Theta  %*% Zt[i-1,]
}
par(mar=c(4,4,2,2))
plot(Xt[,1],type='l',ylim=range(Xt))
lines(Xt[,2],col='red')
@
We can check the ARIMA conclusion by inspecting the ACFs
<<Ex3>>=
cbind(acf(Xt[,1],plot=F,lag.max=10)$acf,acf(Xt[,2],plot=F,lag.max=10)$acf)
cbind(acf(diff(Xt[,1]),plot=F,lag.max=10)$acf,acf(diff(Xt[,2]),plot=F,lag.max=10)$acf)
@ 



However, let
\[
\bL = \left[ \begin{array}{rr} 1 & -2 \\ 0.5 & 1 \end{array} \right]
\]
and let $\bY_t = \bL \bX_t$ and $\bU_t = \bL \bZ_t$
Then
\begin{eqnarray*}
\bY_t = \bL \bX_t & = & \bL \bPhi \bX_{t-1} + \bL \bZ_t + \bL \bTh \bZ_{t-1}\\[6pt]
& = & \bL \bPhi \bL^{-1} (\bL \bX_{t-1}) + \bL \bZ_t + \bL \bTh \bL^{-1} (\bL  \bZ_{t-1})\\[6pt]
& = & \bL \bPhi \bL^{-1} \bY_{t-1} + \bU_t + \bL \bTh \bL^{-1} \bU_{t-1}
\end{eqnarray*}
or equivalently
\[
\bY_t - \bPhi^\star \bY_{t-1} = \bU_t + \bTh^\star \bU_{t-1}
\]

In this equation
\[
\bPhi^\star = \left[ \begin{array}{cc} 1 & 0 \\ 0 & 0 \end{array} \right]
\qquad \qquad
\bTh^\star = \left[ \begin{array}{cc} -0.4 & 0 \\ 0 & 0 \end{array} \right]
\]
so that
\begin{eqnarray*}
Y_{1,t} - Y_{1,t-1} & = & U_{1,t} - 0.4 U_{1,t-1}\\[6pt]
Y_{2,t} & = & U_{2,t}
\end{eqnarray*}
Hence $\{Y_{1t}\} \sim \text{ARIMA}(0,1,1)$ and $\{Y_{2t}\}$ is a white noise process, and there is no dependence between the two series.

<<Ex4,fig.height=5.0,fig.width=7>>=
L<-matrix(c(1,-2,0.5,1),2,2,byrow=T)
Phistar<-L %*% Phi %*% solve(L)
Thetastar<- L %*% Theta %*% solve(L)
Sigmastar<-L %*% Sigma %*% solve(L)
Phistar
Thetastar
Sigmastar

Yt<-t(L %*% t(Xt))
acf.mat<-cbind(acf(Yt[,1],plot=F,lag.max=10)$acf,acf(Yt[,2],plot=F,lag.max=10)$acf,
            acf(diff(Yt[,1]),plot=F,lag.max=10)$acf)
colnames(acf.mat)<-c('Y1','Y2','diff(Y1)')
rownames(acf.mat)<-paste('h=',0:10)
acf.mat
par(mar=c(4,4,3,2))
acf(Yt)

@

\vspace{0.1 in} Hence, the unit roots in $\{X_{1,t} \}$ and $\{X_{2,t} \}$ result from the unit root in $\{Y_{1,t} \}$ and a subsequent linear transformation. This phenomenon is known as \emph{cointegration}.

\pagebreak

\textbf{Example:} UK Consumption Data: In the UK consumption data, the three series are:

\begin{itemize}

\item Consumption, $C_t$

\item Income, $I_t$

\item Wealth, $W_t$

\end{itemize}
in millions of pounds, measured on the log scale.

<<Co0,comment='+', fig.width=7.25, fig.height=5.15, fig.align='center'>>=
library(tseries)
library(urca)
data(Raotbl3)
attach(Raotbl3)
lc <- ts(lc, start=c(1966,4),end=c(1991,2),frequency=4)
li <- ts(li, start=c(1966,4),end=c(1991,2),frequency=4)
lw <- ts(lw, start=c(1966,4),end=c(1991,2),frequency=4)
ukcons <- window(cbind(lc,li,lw), start=c(1967,2),end=c(1991,2))
par(mar=c(4,4,2,2))
ts.plot(cbind(lc,li,lw),ylim=range(10,14),col=c('black','blue','red'))
legend("left", legend = c('Consumption','Income','Wealth'), 
        col = c('black','blue','red'), lty = 1)
@
The three series appear to vary together, but appear to be non-stationary.

\medskip

\textbf{Unit Root Testing:}  The following results show that there seems to be a cointegration relationship between the three series.  We fit the models
\begin{align*}
C_t &= \alpha_0 + \alpha_1 I_t + \alpha_2 W_t + Z_t^\ast\\[6pt]
I_t &= \alpha_0 + \alpha_1 C_t + \alpha_2 W_t + Z_t^\ast\\[6pt]
W_t &= \alpha_0 + \alpha_1 I_t + \alpha_2 C_t + Z_t^\ast
\end{align*}
with the cointegrating vector of $\alpha$s being estimated by OLS.  If the residual series are white noise then there is evidence of cointegration as, for example
\[
C_t - \alpha_0 - \alpha_1 I_t - \alpha_2 W_t
\]
is stationary. The three \texttt{lm} fits below attempt to find the cointegrating vectors.

<<Co1,comment='+', fig.width=7.25, fig.height=5.25, fig.align='center'>>=
lc.eq <- lm(lc~li + lw, data=ukcons)
li.eq <- lm(li~lc + lw, data=ukcons)
lw.eq <- lm(lw~li + lc, data=ukcons)
round(summary(lc.eq)$coef,6)
round(summary(li.eq)$coef,6)
round(summary(lw.eq)$coef,6)
@

The residual series are computed by
<<Co2,comment='+', fig.width=7.25, fig.height=5.75, fig.align='center'>>=
error.lc <- ts(resid(lc.eq), start=c(1967,2),end=c(1991,2), frequency=4)
error.li <- ts(resid(li.eq), start=c(1967,2),end=c(1991,2), frequency=4)
error.lw <- ts(resid(lw.eq), start=c(1967,2),end=c(1991,2), frequency=4)
@

The code below looks for unit roots amongst the residual series:
<<Co3,comment='+', fig.width=7.25, fig.height=8.25, fig.align='center'>>=
ci.lc <- ur.df(error.lc, lags=1, type='none')
ci.li <- ur.df(error.li, lags=1, type='none')
ci.lw <- ur.df(error.lw, lags=1, type='none')

summary(ci.lc)@cval; summary(ci.lc)@teststat
summary(ci.li)@cval; summary(ci.li)@teststat
summary(ci.lw)@cval; summary(ci.lw)@teststat
@
\bigskip
The conclusion is that 

\begin{itemize}[itemsep=0.1in]

\item there appears to be strong evidence for stationarity in the residuals from the regression of $C_t$ on $(I_t,W_t)$ (the test statistic is -4.1414, where the 1\% critical value of -2.6); 

\item there appears to be strong evidence for stationarity in the residuals from the regression of $I_t$ on $(C_t,W_t)$ (the test statistic is -4.059);

\item evidence from the third regression (test statistic -2.7137) is less compelling, but still significant.

\end{itemize}

\bigskip


\textbf{Error Correction Models:}  The three error correction models (ECMs)
<<Co4,comment='+', fig.width=7.25, fig.height=5.25, fig.align='center'>>=
#Construct the error correction models
#Columns contain differenced, and lag-1 differenced data
ukcons2 <- ts(embed(diff(ukcons), dim=2), start=c(1967,4), freq=4)
colnames(ukcons2) <- c('lc.d','li.d','lw.d','lc.d1','li.d1','lw.d1')

#Construct the ECM Error terms (lagged residuals from original regression)
error.ecm1 <- window(lag(error.lc, k=-1),start=c(1967,4), end=c(1991, 2))
error.ecm2 <- window(lag(error.li, k=-1),start=c(1967,4), end=c(1991, 2))
error.ecm3 <- window(lag(error.lw, k=-1),start=c(1967,4), end=c(1991, 2))
ecm.eq1 <- lm(lc.d ~ error.ecm1 + lc.d1 + li.d1 + lw.d1,data=ukcons2)
ecm.eq2 <- lm(li.d ~ error.ecm2 + lc.d1 + li.d1 + lw.d1,data=ukcons2)
ecm.eq3 <- lm(lw.d ~ error.ecm3 + lc.d1 + li.d1 + lw.d1,data=ukcons2)
@
examine the equilibrium relationships between differenced versions of the variables. The three models fitted are
\begin{eqnarray*}
(1-B) C_t &=& \psi_0 + \gamma_1 \zhat_{1,t-1}^\ast + \psi_1 (1-B)C_{t-1} + \psi_2 (1-B)I_{t-1} + \psi_3 (1-B) W_{t-1} + Z_t\\[6pt]
(1-B) I_t &=& \psi_0 + \gamma_2 \zhat_{2,t-1}^\ast + \psi_1 (1-B)C_{t-1} + \psi_2 (1-B)I_{t-1} + \psi_3 (1-B) W_{t-1} + Z_t\\[6pt]
(1-B) W_t &=& \psi_0 + \gamma_3 \zhat_{3,t-1}^\ast + \psi_1 (1-B)C_{t-1} + \psi_2 (1-B)I_{t-1} + \psi_3 (1-B) W_{t-1} + Z_t
\end{eqnarray*}
where $\zhat_{i,t-1}^\ast$ for $i=1,2,3$ are the residuals from the earlier fitted cointegration models.

\bigskip

Estimates of the parameters $\gamma_1,\gamma_2,\gamma_3$ will give information about the ``correction" of the system; for example $\gamma_1$ measures the magnitude of the effect of the (lagged) residual from the series
\[
C_t - (\alpha_0 + \alpha_1 I_t + \alpha_2 W_t) .
\]
If this lagged residual $\zhat_{1,t-1}^\ast$ is positive, then $\gamma_1 < 0$ indicates that the expected \textbf{change} in $C_t$ over the next time interval (as measured by the left-hand side of the ECM, $(1-B)C_t$, will also be negative as $\gamma_1 \zhat_{1,t-1}^\ast < 0$, that is, $C_t$ is ``corrected" by the long-run stability of the system.

\begin{enumerate}[label=(\roman*)]
\item From the first ECM, we obtain the following result:
<<Co5,comment='+', fig.width=7.25, fig.height=5.25, fig.align='center'>>=
summary(ecm.eq1)$coef
@
\bigskip
The coefficient of \texttt{error.ecm1} is not significant, so the error correction from the lagged cointegration relationship is negligible in this case.  Thus there appears to be no ``causal" structure, that is, there is little additional predictive power in understanding the cointegration relationships between the variables.

\bigskip

\item For the second ECM, we have the following results
<<Co6,comment='+', fig.width=7.25, fig.height=5.25, fig.align='center'>>=
summary(ecm.eq2)$coef
@
Here, the coefficient of \texttt{error.ecm2} is significantly different from zero, thus indicating a causal relationship; the estimate is negative as required, indicating a ``correction" mechanism is present. 


\item For the third ECM, there appears to be no causal relationship as the coefficient of \texttt{error.ecm3} is not significantly different from zero.
<<Co7,comment='+', fig.width=7.25, fig.height=5.25, fig.align='center'>>=
summary(ecm.eq3)$coef
@
These findings are confirmed by looking at the correlation between the differenced responses (that is, $(1-B)C_t, (1-B)I_t$, and $(1-B)W_t$) and the cointegration residuals $\zhat_{i,t-1}^\ast$:
<<Co8,comment='+', fig.width=7.25, fig.height=5.25, fig.align='center'>>=
cor(as.numeric(ukcons2[,1]),error.ecm1)
cor(as.numeric(ukcons2[,2]),error.ecm2)
cor(as.numeric(ukcons2[,3]),error.ecm3)
@
\end{enumerate}

\end{document}