\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 in a simple model} \begin{document} <>= library(knitr) # global chunk options opts_chunk$set(cache=TRUE, autodep=TRUE) options(scipen=999) options(repos=c(CRAN="https://cloud.r-project.org/")) myrnd <- function(x,d) {formatC(x, format="f", digits=d)} @ \maketitle Suppose \begin{align*} X_t & = X_{t-1} + Z_{1,t} \\[3pt] Y_t & = 0.6 X_t + Z_{2,t}. \end{align*} Clearly $\{X_t\}$ is non-stationary; also, \[ Y_t - Y_{t-1} = 0.6 (X_t - X_{t-1}) + Z_{2,t} - Z_{2,t-1} = 0.6 Z_{1,t} + Z_{2,t} - Z_{2,t-1} \] so $\{Y_t\}$ is non-stationary, being stationary only after differencing. The model encapsulates a long-run relationship between the two non-stationary series, in that $Y_t - 0.6 X_t$ is stationary. Regressing $Y_t$ on $X_t$ allows the cointegrating relationship to be discovered. \medskip Our task is to predict \textit{changes} in $Y_t$ on the basis of \textit{past} values of the series, that is, data up to and including time point $t-1$. We could try to model changes in $Y_t$ via a linear model based on the stationary differenced series \[ (Y_t - Y_{t-1}) = \psi_0 + \psi_1 (X_{t-1} - X_{t-2}) + \psi_2 (Y_{t-1}-Y_{t-2}) + Z_t \] or \[ Y_t^\ast = \psi_0 + \psi_1 X_{t-1}^\ast + \psi_2 Y_{t-1}^\ast + Z_t \] say. However, this model does not exploit the long-run relationship between $Y_t$ and $X_t$, so we add in a third term to make the \textit{error correction model} \[ Y_t^\ast = \psi_0 + \psi_1 X_{t-1}^\ast + \psi_2 Y_{t-1}^\ast + \gamma (Y_{t-1} - 0.6 X_{t-1}) + Z_t \] that reflects the fact that if indeed there is a long-run relationship between $\{Y_t\}$ and $\{X_t\}$, then a positive discrepancy -- occurring by chance -- at time $t-1$, so that $Y_{t-1} > 0.6 X_{t-1}$, will have to be `corrected' at time $t$. Consequently, the difference $Y_t - Y_{t-1}$ should be smaller (or more negative) than it would be if the equilibrium was observed at time $t-1$ (ie $Y_{t-1} \approx 0.6 X_{t-1}$). \medskip Note that in fact any historical error correction term, that is \[ (Y_{t-j} - \alpha_0 - \alpha_1 X_{t-j}) \qquad j \geq 1 \] could be used in this model to act as an error corrector; the calculation is based on the assumption that the sequence of error correction terms is stationary, and at time $t$ all the historical terms are available. \bigskip \textbf{Implementation} <>= set.seed(1231456) n<-1000 z1<-rnorm(n) z2<-rnorm(n) x<-cumsum(z1); y<-0.6*x+z2 @ We first fit the model \begin{equation}\label{eq:Fit0} Y_t^\ast = \psi_0 + \psi_1 X_{t-1}^\ast + \psi_2 Y_{t-1}^\ast + Z_t \end{equation} <>= ystar<-diff(y)[-1] ystar.1<-diff(y)[-(n-1)] xstar.1<-diff(x)[-(n-1)] fit0<-lm(ystar~xstar.1+ystar.1) summary(fit0$coef) #R-squared summary(fit0)$r.squared @ This shows that the differenced series are useful in forecasting the change in $Y_t$. Now we try the error correction model assuming the known relationship between $\{X_t\}$ and $\{Y_t\}$: \begin{equation}\label{eq:Fit1} Y_t^\ast = \psi_0 + \psi_1 X_{t-1}^\ast + \psi_2 Y_{t-1}^\ast + \gamma (Y_{t-1} - 0.6 X_{t-1}) + Z_t \end{equation} where $(Y_{t-1} - 0.6 X_{t-1}) = Z_{2,t}$. <>= error<-z2; error.lagged<-error[-c(1,n)] xstar<-diff(x);ystar<-diff(y) diff.dat<-data.frame(embed(cbind(xstar,ystar),2)) colnames(diff.dat)<-c("xstar","ystar","xstar.1","ystar.1") ecm.reg0<-lm(ystar~xstar.1+ystar.1+error.lagged,data=diff.dat) round(summary(ecm.reg0)$coef,6) #R-squared summary(ecm.reg0)$r.squared @ This analysis clearly shows that the error correction term is influential, with \[ \widehat \gamma = \Sexpr{myrnd(summary(ecm.reg0)$coef[4,1],5)} \] with standard error \Sexpr{myrnd(summary(ecm.reg0)$coef[4,2],5)}. \bigskip We now consider estimating the model \[ Y_t = \alpha_0 + \alpha_1 X_t + Z_{t} \] to estimate the cointegrating vector and repeating the analysis <>= lr.reg<-lm(y~x) summary(lr.reg)$coef @ We obtain the residuals $\{\widehat{z}_{t}^\star\}$ from the fit, and proceed with the model \begin{equation}\label{eq:Fit2} Y_t^\ast = \psi_0 + \psi_1 X_{t-1}^\ast + \psi_2 Y_{t-1}^\ast + \gamma z_{t-1}^\ast + Z_t \end{equation} <>= error<-residuals(lr.reg) error.lagged<-error[-c(1,n)] xstar<-diff(x);ystar<-diff(y) diff.dat<-data.frame(embed(cbind(xstar,ystar),2)) colnames(diff.dat)<-c("xstar","ystar","xstar.1","ystar.1") ecm.reg1<-lm(ystar~xstar.1+ystar.1+error.lagged,data=diff.dat) round(summary(ecm.reg1)$coef,6) #R-squared summary(ecm.reg1)$r.squared @ Note that the R-squared quantities for the models in equations \eqref{eq:Fit1} and \eqref{eq:Fit2} are higher than that for \eqref{eq:Fit0}, indicating the predictive advantage in adding the error correction term. This analysis clearly shows that the error correction term is influential, with \[ \widehat \gamma = \Sexpr{myrnd(summary(ecm.reg1)$coef[4,1],5)} \] with standard error \Sexpr{myrnd(summary(ecm.reg1)$coef[4,2],5)}. \medskip Finally, for illustration, we demonstrate that any historical error correction term may be included, for example, $z_{t-2}^\ast$ rather than $z_{t-1}^\ast$, using the model \begin{equation}\label{eq:Fit3} Y_t^\ast = \psi_0 + \psi_1 X_{t-1}^\ast + \psi_2 Y_{t-1}^\ast + \gamma z_{t-2}^\ast + Z_t \end{equation} <>= error.lagged<-error[-c(n-1,n)] xstar<-diff(x);ystar<-diff(y) diff.dat<-data.frame(embed(cbind(xstar,ystar),2)) colnames(diff.dat)<-c("xstar","ystar","xstar.1","ystar.1") ecm.reg2<-lm(ystar ~xstar.1+ystar.1+error.lagged,data=diff.dat) round(summary(ecm.reg2)$coef,6) #R-squared summary(ecm.reg2)$r.squared @ This produces an identical fit for $\gamma$, although the estimates for the other regression coefficients change. \bigskip \textbf{Theoretical analysis of the Error Correction Model:} Considering the model in \eqref{eq:Fit1}, we have by standard theory that the least-squares estimators converge to \[ \E[\bX^\top \bX]^{-1} \E[\bX^\top Y] \] where $Y \equiv Y_t^\ast$, and \[ \bX = [1, X_{t-1}^\ast, Y_{t-1}^\ast]. \] is formed from the stationary series. Now \[ \bX^\top \bX = \begin{bmatrix} 1 & X_{t-1}^\ast & Y_{t-1}^\ast\\[6pt] X_{t-1}^\ast & (X_{t-1}^\ast)^2 & X_{t-1}^\ast Y_{t-1}^\ast\\[6pt] Y_{t-1}^\ast & X_{t-1}^\ast Y_{t-1}^\ast & (Y_{t-1}^\ast)^2 \end{bmatrix} \] and we have immediately that $\E[X_{t-1}^\ast] = \E[Y_{t-1}^\ast] = 0$. Also, if $\phi = 0.6$, we have \begin{align*} \E[(X_{t-1}^\ast)^2] & = \E[(X_{t-1}-X_{t-2})^2] = \E[Z_{1,t-1}^2] = \sigma_1^2 \\[6pt] \E[X_{t-1}^\ast Y_{t-1}^\ast ] & = \E[(X_{t-1}-X_{t-2})(\phi (X_{t-1}-X_{t-2})+ Z_{2,t-1}-Z_{2,t-1})] = \phi \E[Z_{1,t-1}^2] = \phi \sigma_1^2\\[6pt] \E[(Y_{t-1}^\ast)^2] & = \E[(\phi (X_{t-1}-X_{t-2})+ Z_{2,t-1}-Z_{2,t-1})^2] = \E[(\phi Z_{1,t-1} + Z_{2,t-1}-Z_{2,t-2})^2] = \phi^2 \sigma_1^2 + 2 \sigma_2^2. \end{align*} Thus \[ \E[\bX^\top \bX] = \begin{bmatrix} 1 & 0 & 0\\[9pt] 0 & \sigma_1^2 & \phi \sigma_1^2\\[9pt] 0 & \phi \sigma_1^2 & \phi^2 \sigma_1^2 + 2 \sigma_2^2. \end{bmatrix} \qquad \therefore \qquad \E[\bX^\top \bX]^{-1} = \begin{bmatrix} 1 & 0 & 0\\[9pt] 0 & \dfrac{\phi^2}{2\sigma_2^2} + \dfrac{1}{2\sigma_1^2}& -\dfrac{\phi}{2\sigma_2^2}\\[9pt] 0 & -\dfrac{\phi}{2\sigma_2^2} & \dfrac{1}{2\sigma_2^2} \end{bmatrix} \] Also \[ \E[\bX^\top Y] = \begin{bmatrix} Y_{t}^\ast \\[3pt] Y_{t}^\ast X_{t-1}^\ast\\[3pt] Y_{t}^\ast Y_{t-1}^\ast \end{bmatrix}= \begin{bmatrix} 0 \\[3pt] 0\\[3pt] -\sigma_2^2 \end{bmatrix} \] as \[ \E[Y_{t}^\ast Y_{t-1}^\ast] = \E[(\phi Z_{1,t} + Z_{2,t}-Z_{2,t-1})((\phi Z_{1,t-1} + Z_{2,t-1}-Z_{2,t-2}))] = -\E[Z_{2,t-1}^2] \] so \[ \E[\bX^\top \bX]^{-1} \E[\bX^\top Y] = \begin{bmatrix} 0 \\[9pt] \dfrac{\phi}{2}\\[9pt] -\dfrac{1}{2} \end{bmatrix} \] \bigskip For the model in \eqref{eq:Fit2}, we have \[ \bX = [1, X_{t-1}^\ast, Y_{t-1}^\ast, Z_{2,t-1}]. \] Now we have \[ \bX^\top \bX = \begin{bmatrix} 1 & X_{t-1}^\ast & Y_{t-1}^\ast & Z_{2,t-1}\\[6pt] X_{t-1}^\ast & (X_{t-1}^\ast)^2 & X_{t-1}^\ast Y_{t-1}^\ast & X_{t-1}^\ast Z_{2,t-1}\\[6pt] Y_{t-1}^\ast & X_{t-1}^\ast Y_{t-1}^\ast & (Y_{t-1}^\ast)^2 & Y_{t-1}^\ast Z_{2,t-1}\\[6pt] Z_{2,t-1} & X_{t-1}^\ast Z_{2,t-1} & Y_{t-1}^\ast Z_{2,t-1} & (Z_{2,t-1})^2 \end{bmatrix} \] where \begin{align*} \E[X_{t-1}^\ast Z_{2,t-1}] & = 0\\[6pt] \E[Y_{t-1}^\ast Z_{2,t-1}] & = \E[(\phi (X_{t-1}-X_{t-2})+ Z_{2,t-1}-Z_{2,t-2}) Z_{2,t-1}] = \E[(\phi Z_{1,t-1} + Z_{2,t-1}-Z_{2,t-2}) Z_{2,t-1}] = \sigma_2^2 \end{align*} Thus \[ \E[\bX^\top \bX] = \begin{bmatrix} 1 & 0 & 0 & 0 \\[9pt] 0 & \sigma_1^2 & \phi \sigma_1^2 & 0 \\[9pt] 0 & \phi \sigma_1^2 & \phi^2 \sigma_1^2 + 2 \sigma_2^2 & \sigma_2^2\\[6pt] 0 & 0 & \sigma_2^2 & \sigma_2^2 \end{bmatrix} \qquad \therefore \qquad \E[\bX^\top \bX]^{-1} = \begin{bmatrix} 1 & 0 & 0 & 0 \\[9pt] 0 & \dfrac{\phi^2}{\sigma_2^2} & -\dfrac{\phi}{\sigma_2^2} & \dfrac{\phi}{\sigma_2^2}\\[9pt] 0 & -\dfrac{\phi}{\sigma_2^2} & \dfrac{1}{\sigma_2^2} & - \dfrac{1}{\sigma_2^2}\\[9pt] 0 & \dfrac{\phi}{\sigma_2^2} & - \dfrac{1}{\sigma_2^2} & \dfrac{2}{\sigma_2^2} \end{bmatrix} \] Also \[ \E[\bX^\top Y] = \begin{bmatrix} Y_{t}^\ast \\[3pt] Y_{t}^\ast X_{t-1}^\ast\\[3pt] Y_{t}^\ast Y_{t-1}^\ast \\[3pt] Y_{t}^\ast Z_{2,t-1} \end{bmatrix}= \begin{bmatrix} 0 \\[3pt] 0\\[3pt] -\sigma_2^2\\[3pt] -\sigma_2^2 \end{bmatrix} \] as \[ \E[Y_{t}^\ast Z_{2,t-1}] = \E[(\phi Z_{1,t} + Z_{2,t}-Z_{2,t-1}) Z_{2,t-1}] = -\E[Z_{2,t-1}^2]. \] Therefore \[ \E[\bX^\top \bX]^{-1} \E[\bX^\top Y] = \begin{bmatrix} 0 \\[6pt] 0\\[6pt] 0\\[6pt] -1 \end{bmatrix} \] \end{document} \medskip Finally, for the model in \eqref{eq:Fit3}, we have \[ \bX = [1, X_{t-1}^\ast, Y_{t-1}^\ast, Z_{2,t-2}]. \] Now we have \[ \bX^\top \bX = \begin{bmatrix} 1 & X_{t-1}^\ast & Y_{t-1}^\ast & Z_{2,t-2}\\[6pt] X_{t-1}^\ast & (X_{t-1}^\ast)^2 & X_{t-1}^\ast Y_{t-1}^\ast & X_{t-1}^\ast Z_{2,t-2}\\[6pt] Y_{t-1}^\ast & X_{t-1}^\ast Y_{t-1}^\ast & (Y_{t-1}^\ast)^2 & Y_{t-1}^\ast Z_{2,t-2}\\[6pt] Z_{2,t-2} & X_{t-1}^\ast Z_{2,t-2} & Y_{t-1}^\ast Z_{2,t-2} & (Z_{2,t-2})^2 \end{bmatrix} \] where \begin{align*} \E[X_{t-1}^\ast Z_{2,t-2}] & = 0\\[6pt] \E[Y_{t-1}^\ast Z_{2,t-2}] & = \E[(\phi (X_{t-1}-X_{t-2})+ Z_{2,t-1}-Z_{2,t-2}) Z_{2,t-2}] = \E[(\phi Z_{1,t-1} + Z_{2,t-1}-Z_{2,t-2}) Z_{2,t-2}] = -\sigma_2^2 \end{align*} Thus \[ \E[\bX^\top \bX] = \begin{bmatrix} 1 & 0 & 0 & 0 \\[9pt] 0 & \sigma_1^2 & \phi \sigma_1^2 & 0 \\[9pt] 0 & \phi \sigma_1^2 & \phi^2 \sigma_1^2 + 2 \sigma_2^2 & -\sigma_2^2\\[6pt] 0 & 0 & -\sigma_2^2 & \sigma_2^2 \end{bmatrix} \qquad \therefore \qquad \E[\bX^\top \bX]^{-1} = \begin{bmatrix} 1 & 0 & 0 & 0 \\[9pt] 0 & \dfrac{\phi^2}{\sigma_2^2} & -\dfrac{\phi}{\sigma_2^2} & \dfrac{\phi}{\sigma_2^2}\\[9pt] 0 & -\dfrac{\phi}{\sigma_2^2} & \dfrac{1}{\sigma_2^2} & - \dfrac{1}{\sigma_2^2}\\[9pt] 0 & \dfrac{\phi}{\sigma_2^2} & - \dfrac{1}{\sigma_2^2} & \dfrac{2}{\sigma_2^2} \end{bmatrix} \] Also \[ \E[\bX^\top Y] = \begin{bmatrix} Y_{t}^\ast \\[3pt] Y_{t}^\ast X_{t-1}^\ast\\[3pt] Y_{t}^\ast Y_{t-1}^\ast \\[3pt] Y_{t}^\ast Z_{2,t-1} \end{bmatrix}= \begin{bmatrix} 0 \\[3pt] 0\\[3pt] -\sigma_2^2\\[3pt] 0 \end{bmatrix} \] as \[ \E[Y_{t}^\ast Z_{2,t-1}] = \E[(\phi Z_{1,t} + Z_{2,t}-Z_{2,t-1}) Z_{2,t-2}] = 0. \] Therefore \[ \E[\bX^\top \bX]^{-1} \E[\bX^\top Y] = \begin{bmatrix} 0 \\[6pt] 0\\[6pt] 0\\[6pt] -1 \end{bmatrix} \]