# CointExamples.ssc script file for examples used in Cointegration chapter # # author: Eric Zivot # created: December 26, 2001 # revised: March 23, 2005 # # Notes: ################################################################### # Spurious regression set.seed(458) e1 = rnorm(250) e2 = rnorm(250) y1 = cumsum(e1) y2 = cumsum(e2) tsplot(y1,y2,lty=c(1,3)) legend(0,15,c("y1","y2"),lty=c(1,3)) summary(OLS(y1~y2)) summary(OLS(diff(y1)~diff(y2))) # Example 1: simulate cointegrated processes # Phillips' triangular representation # 1 cointegrating vector, 1 common trend # x.rw and y.ur are cointegrated with cointegrating vector (1,-1) # and ar1 errors. x.rw is the common trend. # set.seed(432) e = rmvnorm(250,mean=rep(0,2),sd=c(0.5,0.5)) u.ar1 = arima.sim(model=list(ar=0.75),innov=e[,1]) y2 = cumsum(e[,2]) y1 = y2 + u.ar1 par(mfrow=c(2,1)) tsplot(y1,y2,main="Simulated bivariate cointegrated system", sub="1 cointegrating vector, 1 common trend",lty=c(1,3)) legend(0,7,legend=c("y1","y2"),lty=c(1,3)) tsplot(u.ar1,main="Cointegrating residual") par(mfrow=c(1,1)) # create dataframe for later use coint.dat1 = cbind(y1,y2) dimnames(coint.dat1) = list(NULL,c("y1","y2")) coint.dat1 = data.frame(coint.dat1) # Example 2: simulated cointegrated process # Phillips' triangular representation # 1 cointegrating vector, 2 common trends # cointegrating vector is beta' = (1,-0.5,-0.5) # y.ur is cointegrated with x.rw and z.rw # x.rw is one common trend # z.rw is the second common trend # set.seed(573) e = rmvnorm(250,mean=rep(0,3),sd=c(0.5,0.5,0.5)) u1.ar1 = arima.sim(model=list(ar=0.75),innov=e[,1]) y2 = cumsum(e[,2]) y3 = cumsum(e[,3]) y1 = 0.5*y2 + 0.5*y3 + u1.ar1; par(mfrow=c(2,1)) tsplot(y1,y2,y3,main="Simulated trivariate cointegrated system", sub="1 cointegrating vector, 2 common trends",lty=c(1,3,4)) legend(0,12,legend=c("y1","y2","y3"),lty=c(1,3,4)) tsplot(u.ar1,main="Cointegrating residual") par(mfrow=c(1,1)) # create dataframe for later use coint.dat2 = cbind(y1,y2,y3) dimnames(coint.dat2) = list(NULL,c("y1","y2","y3")) coint.dat2 = data.frame(coint.dat2) # Example 3: simulated cointegrated process # Phillips' triangular representation # 2 cointegrating vector, 1 common trend # x.rw and y1.ur are cointegrated with cointegrating vector (1,0,-1) # and ar1 errors. # x.rw and y2.ur and cointegrated with cointegrating vector (0, 1, -1) # x.rw is the common trend. # set.seed(573) e = rmvnorm(250,mean=rep(0,3),sd=c(0.5,0.5,0.5)) u.ar1 = arima.sim(model=list(ar=0.75),innov=e[,1]) v.ar1 = arima.sim(model=list(ar=0.75),innov=e[,2]) y3 = cumsum(e[,3]) y1 = y3 + u.ar1 y2 = y3 + v.ar1 par(mfrow=c(2,1)) tsplot(y1,y2,y3,main="Simulated trivariate cointegrated system", sub="2 cointegrating vectors, 1 common trend",lty=c(1,3,4)) legend(0,10,legend=c("y1","y2","y3"),lty=c(1,3,4)) tsplot(u.ar1,v.ar1,main="Cointegrated residuals", lty=c(1,3)) legend(0,-1,legend=c("u","v"),lty=c(1,3)) par(mfrow=c(1,1)) # create dataframe for later use coint.dat3 = cbind(y1,y2,y3) dimnames(coint.dat3) = list(NULL,c("y1","y2","y3")) coint.dat3 = data.frame(coint.dat3) # # # test for cointegration between spot and forward exchange rates # using a known cointegrating vector # uscn.s = lexrates.dat[,"USCNS"] uscn.s@title = "Log of US/CA spot exchange rate" uscn.f = lexrates.dat[,"USCNF"] uscn.f@title = "Log of US/CA 30-day forward exchange rate" u = uscn.s - uscn.f colIds(u) = "USCNID" u@title = "US/CA 30-day interest rate differential" par(mfrow=c(2,1)) plot(uscn.s,uscn.f,main="Log US/CA exchange rate data", plot.args=list(lty=c(1,3))) legend(0.8,0,legend=c("spot","forward"),lty=c(1,3)) plot(u) par(mfrow=c(1,1)) unitroot(u,trend="c",method="adf",lags=11) # # compute Horvath-Watson Wald statistic # d.uscn.s = diff(uscn.s) colIds(d.uscn.s) = "DUSCNS" d.uscn.f = diff(uscn.f) colIds(d.uscn.f) = "DUSCNF" u.lag = tslag(u,trim=T) uscn.dat = seriesMerge(d.uscn.s,d.uscn.f,u.lag) colIds(uscn.dat) formula.list = list(DUSCNS~USCNID.lag1, DUSCNF~USCNID.lag1) sur.fit = SUR(formula.list,data=uscn.dat) vecb = as.vector(coef(sur.fit)) bigR = matrix(c(0,1,0,0, 0,0,0,1),2,4,byrow=T) avar = bigR%*%vcov(sur.fit)%*%t(bigR) wald = t(bigR%*%vecb)%*%solve(avar)%*%(bigR%*%vecb) wald # compare pcoint with punitroot args(pcoint) args(qcoint) # PO quantiles adjusted for constant for ADF t-stat with n-1=3 qcoint(c(0.1,0.05,0.01),n.sample=100,n.series=4, trend="c",statistic="t") qunitroot(c(0.1,0.05,0.01),n.sample=100, trend="c",statistic="t") # # test for cointegration between log stock prices and log dividends # colIds(shiller.annual) ln.p = log(shiller.annual[,"real.price"]) colIds(ln.p) = "ln.p" ln.d = log(shiller.annual[,"real.dividend"]) colIds(ln.d) = "ln.d" ln.dpratio = ln.d - ln.p colIds(ln.dpratio) = "ln.dpratio" stockdiv.ts = seriesMerge(ln.p,ln.d,ln.dpratio) stockdiv.ts@title = "Log stock price and dividend data" stockdiv.ts@documentation = c("Natural log of real stock prices and real dividends", "Stock prices are on the Standard & Poor's composite index", "source: Robert Shiller's homepage,", "http://aida.econ.yale.edu/~shiller/data.htm") par(mfrow=c(2,1)) plot(ln.p,ln.d,main="Log stock prices and dividends", plot.args=list(lty=c(1,3))) legend(0,7,legend=c("stock prices","dividends"),lty=c(1,3)) plot(ln.dpratio,main="log dividend price ratio") par(mfrow=c(1,1)) smpl = (positions(ln.dpratio) <= timeDate("1/1/1995")) summary(unitroot(ln.dpratio[smpl,],lags=3)) # # test for cointegration between spot and forward exchange rates # using an estimated cointegrating vector # uscn.ts = seriesMerge(uscn.s,uscn.f) ols.fit = OLS(USCNS~USCNF,data=uscn.ts) ols.fit u.hat = residuals(ols.fit) adf.fit = unitroot(u.hat,trend="nc",method="adf",lags=11) adf.tstat = adf.fit$sval adf.tstat pp.fit = unitroot(u.hat,trend="nc",method="pp") pp.tstat = pp.fit$sval pp.tstat # compare pcoint with punitroot qcoint(c(0.10,0.05,0.01),n.sample=nrow(uscn.s),n.series=2, trend="c",statistic="t") pcoint(adf.tstat,n.sample=nrow(uscn.s),n.series=2, trend="c",statistic="t") pcoint(pp.tstat,n.sample=nrow(uscn.s),n.series=2, trend="c",statistic="t") # # estimate cointegrating vector for exchange rate data using stock and Watson estimator # uscn.df = diff(uscn.f) colIds(uscn.df) = "D.USCNF" uscn.df.lags = tslag(uscn.df,-3:3,trim=T) uscn.ts = seriesMerge(uscn.s,uscn.f,uscn.df.lags) colIds(uscn.ts) dols.fit = OLS(USCNS~USCNF +D.USCNF.lead3+D.USCNF.lead2+D.USCNF.lead1 +D.USCNF.lag0+D.USCNF.lag1+D.USCNF.lag2+D.USCNF.lag3, data=uscn.ts) summary(dols.fit,correction="nw") dols.fit = OLS(tslag(USCNS,-1)~USCNF +D.USCNF.lead3+D.USCNF.lead2+D.USCNF.lead1 +D.USCNF.lag0+D.USCNF.lag1+D.USCNF.lag2+D.USCNF.lag3, data=uscn.ts,na.rm=T) summary(dols.fit,correction="nw") # compare to OLS summary(OLS(tslag(USCNS,-1)~USCNF,data=uscn.ts,na.rm=T)) # # estimate cointegrating vector for stock price and dividend data using # stock and Watson estimator # dln.d = diff(ln.d) colIds(dln.d) = "dln.d" dln.d.lags = tslag(dln.d,-3:3,trim=T) sd.ts = seriesMerge(ln.p,ln.d,dln.d.lags) colIds(sd.ts) dols.fit = OLS(ln.p~ln.d +dln.d.lead1+dln.d.lead2+dln.d.lead3 +dln.d.lag0+dln.d.lag1+dln.d.lag2+dln.d.lag3, data=sd.ts,end=timeDate("1/1/1995")) lrv.hat = asymp.var(residuals(dols.fit),bandwidth=5) varu.hat = sum(residuals(dols.fit)^2)/dols.fit$df.resid ratio = sqrt(lrv.hat/varu.hat) ratio summary(dols.fit) summary(dols.fit,correction="nw") ols.fit = OLS(ln.p~ln.d,data=stockdiv.ts,end=timeDate("1/1/1995")) summary(ols.fit) # # estimate error correction model for exchange rate data # u.hat = uscn.s - 1.004*uscn.f colIds(u.hat) = "U.HAT" uscn.ds = diff(uscn.s) colIds(uscn.ds) = "D.USCNS" uscn.df = diff(uscn.f) colIds(uscn.df) = "D.USCNF" uscn.ts = seriesMerge(uscn.s,uscn.f,uscn.ds,uscn.df,u.hat) ecm.s.fit = OLS(D.USCNS~tslag(U.HAT)+tslag(D.USCNS) +tslag(D.USCNF),data=uscn.ts,na.rm=T) ecm.f.fit = OLS(D.USCNF~tslag(U.HAT)+tslag(D.USCNS)+tslag(D.USCNF), data=uscn.ts,na.rm=T) summary(ecm.s.fit) summary(ecm.f.fit) # # estimate ecm for stock price and earnings data from Shiller # dln.d = diff(ln.d) colIds(dln.d) = "dln.d" dln.p = diff(ln.p) colIds(dln.p) = "dln.p" stockdiv.ts = seriesMerge(dln.p,dln.d,ln.dpratio) ecm.p.fit = OLS(dln.p~tslag(ln.dpratio)+tslag(dln.p)+tslag(dln.d), data=stockdiv.ts,end=timeDate("1/1/1995"),na.rm=T) summary(ecm.p.fit) ecm.d.fit = OLS(dln.d~tslag(ln.dpratio)+tslag(dln.p)+tslag(dln.d), data=stockdiv.ts,end=timeDate("1/1/1995"),na.rm=T) summary(ecm.d.fit) # # simulate cointegrated processes for the 5 trend cases # # # simulate cointegrated VAR(1) model # beta' = (1,-1), alpha' = (-0.5,0) # set.seed(578) nobs = 250 e.var = rmvnorm(nobs,sd=c(1,1),rho=0.5) e.var = t(e.var) beta = c(1,-1) alpha = c(-0.5,0) alpha.perp = c(0,1) mu0 = c(4,0.5) mu1 = c(0.05,0.1) # solve for orthogonal decomposition parameters rho0 = solve(t(alpha)%*%alpha)%*%t(alpha)%*%mu0 rho1 = solve(t(alpha)%*%alpha)%*%t(alpha)%*%mu1 gamma0 = solve(t(alpha.perp)%*%alpha.perp)%*%t(alpha.perp)%*%mu0 gamma1 = solve(t(alpha.perp)%*%alpha.perp)%*%t(alpha.perp)%*%mu1 # check calculations check.mu0 = (mu0 == (alpha%*%rho0 + alpha.perp%*%gamma0)) check.mu1 = (mu1 == (alpha%*%rho1 + alpha.perp%*%gamma1)) check.mu0 check.mu1 # solve for VAR parameters BigPi = as.matrix(alpha)%*%beta pi1 = diag(2) + BigPi # case 1: no constant and trend mu0 = 0 mu1 = 0 y.case1 = matrix(0,2,nobs) for (i in 2:nobs) { y.case1[,i] = mu0+(mu1*i)+pi1%*%y.case1[,i-1]+e.var[,i] } y.case1 = t(y.case1) dimnames(y.case1) = list(NULL,c("y1","y2")) u.case1 = y.case1[,1]-y.case1[,2] # case 2: restricted constant mu0 = as.matrix(alpha%*%rho0) mu1 = 0 y.case2 = matrix(0,2,nobs) for (i in 2:nobs) { y.case2[,i] = mu0+(mu1*i)+pi1%*%y.case2[,i-1]+e.var[,i] } y.case2 = t(y.case2) dimnames(y.case2) = list(NULL,c("y1","y2")) u.case2 = y.case2[,1]-y.case2[,2] # case 3: unrestricted constant mu0 = as.matrix(alpha%*%rho0)+as.matrix(alpha.perp%*%gamma0) mu1 = 0 y.case3 = matrix(0,2,nobs) for (i in 2:nobs) { y.case3[,i] = mu0+(mu1*i)+pi1%*%y.case3[,i-1]+e.var[,i] } y.case3 = t(y.case3) dimnames(y.case3) = list(NULL,c("y1","y2")) u.case3 = y.case3[,1]-y.case3[,2] # question: why doesn't u.case3 have mean 1? # case 4: restricted trend mu0 = as.matrix(alpha%*%rho0)+as.matrix(alpha.perp%*%gamma0) mu1 = as.matrix(alpha%*%rho1) y.case4 = matrix(0,2,nobs) for (i in 2:nobs) { y.case4[,i] = mu0+(mu1*i)+pi1%*%y.case4[,i-1]+e.var[,i] } y.case4 = t(y.case4) dimnames(y.case4) = list(NULL,c("y1","y2")) u.case4 = y.case4[,1]-y.case4[,2] # case 5: unrestricted constant and trend mu0 = as.matrix(alpha%*%rho0)+as.matrix(alpha.perp%*%gamma0) mu1 = as.matrix(alpha%*%rho1)+as.matrix(alpha.perp%*%gamma1) y.case5 = matrix(0,2,nobs) for (i in 2:nobs) { y.case5[,i] = mu0+(mu1*i)+pi1%*%y.case5[,i-1]+e.var[,i] } y.case5 = t(y.case5) dimnames(y.case5) = list(NULL,c("y1","y2")) u.case5 = y.case5[,1]-y.case5[,2] # plot all of the cases smpl = 50:250 par(mfrow=c(3,2)) tsplot(y.case1[smpl,],main="Case 1: No constant") tsplot(y.case2[smpl,],main="Case 2: Restricted constant") tsplot(y.case3[smpl,],main="Case 3: Unrestricted constant") tsplot(y.case4[smpl,],main="Case 4: Restricted trend") tsplot(y.case5[smpl,],main="Case 5: Unrestricted trend") par(mfrow=c(1,1)) par(mfrow=c(3,2)) tsplot(u.case1[smpl],main="Case 1: No constant") tsplot(u.case2[smpl],main="Case 2: Restricted constant") tsplot(u.case3[smpl],main="Case 3: Unrestricted constant ") tsplot(u.case4[smpl],main="Case 4: Restricted trend") tsplot(u.case5[smpl],main="Case 5: Unrestricted trend") par(mfrow=c(1,1)) # # Testing for cointegration when the cointegrating vector is known # # use spot and forward rate data between US and Canada and PP test uscn.fe = lexrates.dat[,"USCNS"] - tslag(lexrates.dat[,"USCNF"]) colIds(uscn.fe) = "USCNFE" uscn.fe@title = "US/CA Excess Return" plot(uscn.fe) pp.test = unitroot(uscn.fe,trend="c",method="pp",na.rm=T) pp.test # use P/E ratio in logs # # Estimate Cointegrating Regression by OLS and DOLS # ols.fit = OLS(USUKS~tslag(USUKF),data=lexrates.dat,na.rm=T) summary(ols.fit) # compute t-stat for testing beta=1 (coef(ols.fit)[2]-1)/sqrt(diag(vcov(ols.fit)))[2] # estimate Stock-Watson DOLS # create data set for example usuks = lexrates.dat[,"USUKS"] usukf = tslag(lexrates.dat[,"USUKF"],1:2) usukf = usukf[,1] colIds(usukf) = "USUKF" dusukf = diff(usukf,trim=F) colIds(dusukf) = "DUSUKF" usuk.dat = seriesMerge(usuks,usukf,dusukf) dols.fit = OLS(USUKS~USUKF+tslag(DUSUKF,-3:3), data=usuk.dat,na.rm=T) summary(dols.fit,correction="nw") # # Johansen cointegration tests # args(coint) # # Exchange rate data # uscn.ts = lexrates.dat[,c("USCNS","USCNF")] # determine lag length using AIC var.fit = VAR(uscn.ts,max.ar=6,criterion="AIC") var.fit$ar.order coint.rc = coint(uscn.ts,trend="rc",lags=1) class(coint.rc) names(coint.rc) coint.rc summary(coint.rc) # # testing linear restrictions # colIds(johansen.danish) seriesPlot(johansen.danish,one.plot=F,type="l") # test if Rb = 0 lies in cointegration space R = c(1, 1, 0, 0, 0) H = perpMat(R) H restr.mod1 = coint(johansen.danish[,c(1,2,4,5)], trend="rc", H=H) print(restr.mod1, restrictions=T) # test b = (1,-1,0,0,0) b = as.matrix(c(1,-1,0,0,0)) restr.mod2 = coint(johansen.danish[,c(1,2,4,5)], trend="rc", b=b) print(restr.mod2, restrictions=T) # un-normalized and normalized cointegrating vectors with r=1 coint.rc$coint.vectors[1,] coint.rc$coint.vectors[1,]/ as.numeric(-coint.rc$coint.vectors[1,1]) # unrestricted constant coint.c = coint(uscn.ts,trend="c",lags=1) vecm.c.fit = VECM(coint.c,unbiased=F) summary(vecm.c.fit) # vecm estimation args(VECM) args(VECM.coint) args(VECM.default) vecm.fit = VECM(coint.rc) class(vecm.fit) inherits(vecm.fit,"VAR") names(vecm.fit) vecm.fit$beta.c vecm.fit$se.beta # note: R^2 values are different from Eviews due to different way they # are calculated for models without a constant vecm.fit summary(vecm.fit) plot(vecm.fit, ask=F) # VECM with prespecified cointegrating vector beta.true = as.matrix(c(1,-1, 0)) dimnames(beta.true) = list(c("USCNS","USCNF","Intercept"), "coint.1") vecm.fruh.fit = VECM(uscn.ts, coint.vec = beta.true, lags = 1, trend = "rc") summary(vecm.fruh.fit) # forecast from vecm vecm.fcst = predict(vecm.fit,n.predict=12) class(vecm.fcst) summary(vecm.fcst) plot(vecm.fcst,xold=diff(uscn.ts),n.old=12) vecm.fit.level = VECM(coint.rc,levels=T) vecm.fcst.level = predict(vecm.fit.level,n.predict=12) summary(vecm.fcst.level) plot(vecm.fcst.level,xold=uscn.ts,n.old=12)