#2 LINEAR MODELS FOR FINANCIAL TIME SERIES
#Chapter 2: Exercise 1,7,8
#Exercise 1
#(a)UnitRoots
library(fUnitRoots)
## Loading required package: timeDate
## Loading required package: timeSeries
## Loading required package: fBasics
da=read.table("m-unrate-4811.txt",header=T)
gdp=log(da[,1])
m1=ar(diff(gdp),method="mle")
m1$order
## [1] 12
adfTest(gdp,lags=2,type=c("c"))
##
## Title:
## Augmented Dickey-Fuller Test
##
## Test Results:
## PARAMETER:
## Lag Order: 2
## STATISTIC:
## Dickey-Fuller: -2.6092
## P VALUE:
## 0.09347
##
## Description:
## Thu Oct 17 23:33:24 2019 by user: TS
#test shows that p-value:0.09347>than significance level 5% which means
#can't reject the null hypothesis. In other words, data has unit root.
#Exercise 1(b)
#2.7.3 Trend-Stationary Time Series
library(fUnitRoots)
da=read.table("m-unrate-4811.txt.",header=T)
unr=log(da[,1])
m2=ar(diff(unr),method='mle') # Based on AIC
m2$order
## [1] 12
adfTest(unr,lags=2,type=("ct"))
##
## Title:
## Augmented Dickey-Fuller Test
##
## Test Results:
## PARAMETER:
## Lag Order: 2
## STATISTIC:
## Dickey-Fuller: -2.8928
## P VALUE:
## 0.2004
##
## Description:
## Thu Oct 17 23:33:25 2019 by user: TS
adfTest(unr,lags=15,type=("ct")) # Based on PACF
##
## Title:
## Augmented Dickey-Fuller Test
##
## Test Results:
## PARAMETER:
## Lag Order: 15
## STATISTIC:
## Dickey-Fuller: -3.202
## P VALUE:
## 0.08758
##
## Description:
## Thu Oct 17 23:33:25 2019 by user: TS
dunr=diff(unr)
tdx=c(1:length(dunr))
m3=arima(dunr,order=c(2,0,0),xreg=tdx)
m3
##
## Call:
## arima(x = dunr, order = c(2, 0, 0), xreg = tdx)
##
## Coefficients:
## ar1 ar2 intercept tdx
## 0.1056 0.2374 0.0018 0
## s.e. 0.0352 0.0353 0.0041 0
##
## sigma^2 estimated as 0.001415: log likelihood = 1425.87, aic = -2841.75
m3$coef
## ar1 ar2 intercept tdx
## 1.055739e-01 2.374440e-01 1.803328e-03 -1.307073e-06
sqrt(diag(m3$var.coef))
## ar1 ar2 intercept tdx
## 3.523565e-02 3.526080e-02 4.130150e-03 3.726187e-05
#exponential smoothing
da=read.table("m-unrate-4811.txt",header=T)
umr=log(da[,1])
length(umr)
## [1] 767
m1=arima(umr,order=c(0,1,1))
m1
##
## Call:
## arima(x = umr, order = c(0, 1, 1))
##
## Coefficients:
## ma1
## 0.0958
## s.e. 0.0305
##
## sigma^2 estimated as 0.001509: log likelihood = 1401.28, aic = -2798.56
Box.test(m1$residuals,lag=10,type='Ljung')
##
## Box-Ljung test
##
## data: m1$residuals
## X-squared = 116.25, df = 10, p-value < 2.2e-16
pp=1-pchisq(14.25,9)
pp
## [1] 0.113706
#forecasting
da=read.table("m-unrate-4811.txt",header=T)
head(da)
## rate
## 1 3.4
## 2 3.8
## 3 4.0
## 4 3.9
## 5 3.5
## 6 3.6
umr=log(da[,1])
m1=arima(umr,order=c(0,0,9)) # unrestricted model
m1
##
## Call:
## arima(x = umr, order = c(0, 0, 9))
##
## Coefficients:
## ma1 ma2 ma3 ma4 ma5 ma6 ma7 ma8
## 1.3266 1.6191 1.7721 1.7477 1.6832 1.4583 1.1231 0.7387
## s.e. 0.0375 0.0628 0.0827 0.0898 0.0794 0.0667 0.0601 0.0498
## ma9 intercept
## 0.3636 1.7092
## s.e. 0.0321 0.0202
##
## sigma^2 estimated as 0.001914: log likelihood = 1309.46, aic = -2596.92
m1=arima(umr,order=c(0,0,9),fixed=c(NA,0,NA,0,0,0,0,0,NA,NA))
m1
##
## Call:
## arima(x = umr, order = c(0, 0, 9), fixed = c(NA, 0, NA, 0, 0, 0, 0, 0, NA, NA))
##
## Coefficients:
## ma1 ma2 ma3 ma4 ma5 ma6 ma7 ma8 ma9 intercept
## 1.1759 0 0.6646 0 0 0 0 0 1.0573 1.7115
## s.e. 0.0752 0 0.0469 0 0 0 0 0 0.0582 0.0115
##
## sigma^2 estimated as 0.006772: log likelihood = 508.32, aic = -1006.64
sqrt(0.005097)
## [1] 0.07139328
Box.test(m1$residuals,lag=12,type='Ljung')# model checking
##
## Box-Ljung test
##
## data: m1$residuals
## X-squared = 4723.9, df = 12, p-value < 2.2e-16
pv=1-pchisq(17.6,9) # compute p-value after adjusting the d.f.
pv
## [1] 0.04010828
m1=arima(umr[1:986],order=c(0,0,9),fixed=c(NA,0,NA,0,0,0,0,0,NA,NA))
m1
##
## Call:
## arima(x = umr[1:986], order = c(0, 0, 9), fixed = c(NA, 0, NA, 0, 0, 0, 0, 0,
## NA, NA))
##
## Coefficients:
## ma1 ma2 ma3 ma4 ma5 ma6 ma7 ma8 ma9 intercept
## 1.1758 0 0.6646 0 0 0 0 0 1.0574 1.7115
## s.e. 0.0752 0 0.0469 0 0 0 0 0 0.0582 0.0115
##
## sigma^2 estimated as 0.006773: log likelihood = 508.32, aic = -1006.64
predict(m1,4) # prediction
## Warning in predict.Arima(m1, 4): MA part of model is not invertible
## $pred
## Time Series:
## Start = 987
## End = 990
## Frequency = 1
## [1] 1.711524 1.711524 1.711524 1.711524
##
## $se
## Time Series:
## Start = 987
## End = 990
## Frequency = 1
## [1] 0.1633997 0.1633997 0.1633997 0.1633997
#Exercise 1(c)
da=read.table("m-unrate-4811.txt",header=T)
head(da)
## rate
## 1 3.4
## 2 3.8
## 3 4.0
## 4 3.9
## 5 3.5
## 6 3.6
umr=log(da[,1])
kumr=ts(umr,frequency=4,start=c(1983,1))
plot(kumr,type='l')
c1=c("2","3","4","1")
c2=c("1","2","3","4")
points(kumr,pch=c1,cex=0.6)
par(mfcol=c(2,2))
kumr=log(da[,1])
dumr=diff(kumr)
skumr=diff(kumr,4)
dkumr=diff(skumr)
acf(kumr,lag=10)
acf(dumr,lag=10)
acf(skumr,lag=10)
acf(dkumr,lag=10)
par(mfcol=c(3,1))
plot(dumr,xlab='month',ylab='diff',type='l')
points(dumr,pch=c1,cex=0.7)
plot(skumr,xlab='month',ylab='sea-diff',type='l')
points(skumr,pch=c2,cex=0.7)
plot(dkumr,xlab='month',ylab='dd',type='l')
points(dkumr,pch=c1,cex=0.7)
# Estimation
m1=arima(kumr,order=c(0,1,1),seasonal=list(order=c(0,1,1),period=4))
m1
##
## Call:
## arima(x = kumr, order = c(0, 1, 1), seasonal = list(order = c(0, 1, 1), period = 4))
##
## Coefficients:
## ma1 sma1
## 0.0955 -1.0000
## s.e. 0.0306 0.0074
##
## sigma^2 estimated as 0.001514: log likelihood = 1382.06, aic = -2758.12
tsdiag(m1,gof=20) # model checking
Box.test(m1$residuals,lag=12,type='Ljung')
##
## Box-Ljung test
##
## data: m1$residuals
## X-squared = 143.2, df = 12, p-value < 2.2e-16
pp=1-pchisq(13.30,10)
pp
## [1] 0.2073788
# Out-of-sample forecasting
kumr=log(da[,1])
length(kumr)
## [1] 767
y=kumr[1:100]
m1=arima(y,order=c(0,1,1),seasonal=list(order=c(0,1,1),period=4))
m1
##
## Call:
## arima(x = y, order = c(0, 1, 1), seasonal = list(order = c(0, 1, 1), period = 4))
##
## Coefficients:
## ma1 sma1
## 0.1522 -0.9996
## s.e. 0.0923 0.2234
##
## sigma^2 estimated as 0.004837: log likelihood = 112.04, aic = -218.08
# Prediction
pm1=predict(m1,10)
names(pm1)
## [1] "pred" "se"
pred=pm1$pred
##Exercise 7
vw=read.table('m-ibm3dx2608.txt',header=T)[,3]
m3=arima(vw,order=c(3,0,0))
m3
##
## Call:
## arima(x = vw, order = c(3, 0, 0))
##
## Coefficients:
## ar1 ar2 ar3 intercept
## 0.1158 -0.0187 -0.1042 0.0089
## s.e. 0.0315 0.0317 0.0317 0.0017
##
## sigma^2 estimated as 0.002875: log likelihood = 1500.86, aic = -2991.73
(1-.1158+.0187+.1042)*mean(vw) # Compute the intercept phi(0).
## [1] 0.008967611
sqrt(m3$sigma2)
## [1] 0.0536189
Box.test(m3$residuals,lag=12,type='Ljung')
##
## Box-Ljung test
##
## data: m3$residuals
## X-squared = 16.352, df = 12, p-value = 0.1756
pv=1-pchisq(16.35,9) # Compute p value using 9 degrees of freedom
pv
## [1] 0.05992276
m3=arima(vw,order=c(3,0,0),fixed=c(NA,0,NA,NA))
## Warning in arima(vw, order = c(3, 0, 0), fixed = c(NA, 0, NA, NA)): some AR
## parameters were fixed: setting transform.pars = FALSE
m3
##
## Call:
## arima(x = vw, order = c(3, 0, 0), fixed = c(NA, 0, NA, NA))
##
## Coefficients:
## ar1 ar2 ar3 intercept
## 0.1136 0 -0.1063 0.0089
## s.e. 0.0313 0 0.0315 0.0017
##
## sigma^2 estimated as 0.002876: log likelihood = 1500.69, aic = -2993.38
(1-.1136+.1063)*.0089
## [1] 0.00883503
sqrt(m3$sigma2)
## [1] 0.05362832
Box.test(m3$residuals,lag=12,type='Ljung')
##
## Box-Ljung test
##
## data: m3$residuals
## X-squared = 16.828, df = 12, p-value = 0.1562
pv=1-pchisq(16.83,10)
pv
## [1] 0.07821131
###############
da=read.table("q-jnj-earns-9211.txt",header=T)
head(da)
## day mon year earns
## 1 30 1 1992 0.11
## 2 23 4 1992 0.18
## 3 21 7 1992 0.18
## 4 20 10 1992 0.17
## 5 1 2 1993 0.12
## 6 29 4 1993 0.20
jnj=log(da[,3])
m1=arima(jnj,order=c(0,0,9)) # unrestricted model
m1
##
## Call:
## arima(x = jnj, order = c(0, 0, 9))
##
## Coefficients:
## ma1 ma2 ma3 ma4 ma5 ma6 ma7 ma8
## 1.9057 2.8842 3.6634 4.9178 4.4922 3.5449 2.5684 1.7858
## s.e. 0.1183 0.2086 0.3142 0.3802 0.4077 0.3567 0.2734 0.1750
## ma9 intercept
## 0.6038 7.6016
## s.e. 0.1004 0.0005
##
## sigma^2 estimated as 2.685e-08: log likelihood = 553.68, aic = -1085.37
m1=arima(jnj,order=c(0,0,9),fixed=c(NA,0,NA,0,0,0,0,0,NA,NA))
m1
##
## Call:
## arima(x = jnj, order = c(0, 0, 9), fixed = c(NA, 0, NA, 0, 0, 0, 0, 0, NA, NA))
##
## Coefficients:
## ma1 ma2 ma3 ma4 ma5 ma6 ma7 ma8 ma9 intercept
## 1.1817 0 0.9295 0 0 0 0 0 1.1196 7.6015
## s.e. 0.2847 0 0.2871 0 0 0 0 0 0.2394 0.0004
##
## sigma^2 estimated as 5.592e-07: log likelihood = 409.57, aic = -809.14
sqrt(0.005097)
## [1] 0.07139328
Box.test(m1$residuals,lag=12,type='Ljung')# model checking
##
## Box-Ljung test
##
## data: m1$residuals
## X-squared = 467.08, df = 12, p-value < 2.2e-16
pv=1-pchisq(17.6,9) # compute p-value after adjusting the d.f.
pv
## [1] 0.04010828
m1=arima(jnj[1:986],order=c(0,0,9),fixed=c(NA,0,NA,0,0,0,0,0,NA,NA))
m1
##
## Call:
## arima(x = jnj[1:986], order = c(0, 0, 9), fixed = c(NA, 0, NA, 0, 0, 0, 0, 0,
## NA, NA))
##
## Coefficients:
## ma1 ma2 ma3 ma4 ma5 ma6 ma7 ma8 ma9 intercept
## 1.1797 0 0.9277 0 0 0 0 0 1.1181 7.6015
## s.e. 0.2836 0 0.2859 0 0 0 0 0 0.2389 0.0004
##
## sigma^2 estimated as 5.606e-07: log likelihood = 409.57, aic = -809.14
predict(m1,10) # prediction
## Warning in predict.Arima(m1, 10): MA part of model is not invertible
## $pred
## Time Series:
## Start = 987
## End = 996
## Frequency = 1
## [1] 7.601529 7.601529 7.601529 7.601529 7.601529 7.601529 7.601529
## [8] 7.601529 7.601529 7.601529
##
## $se
## Time Series:
## Start = 987
## End = 996
## Frequency = 1
## [1] 0.001588704 0.001588704 0.001588704 0.001588704 0.001588704
## [6] 0.001588704 0.001588704 0.001588704 0.001588704 0.001588704
da=read.table("q-jnj-earns-9211.txt",header=T)
head(da)
## day mon year earns
## 1 30 1 1992 0.11
## 2 23 4 1992 0.18
## 3 21 7 1992 0.18
## 4 20 10 1992 0.17
## 5 1 2 1993 0.12
## 6 29 4 1993 0.20
jnj=log(da[,3])
kjnj=ts(jnj,frequency=4,start=c(1983,1))
plot(kjnj,type='l')
c1=c("2","3","4","1")
c2=c("1","2","3","4")
points(kjnj,pch=c1,cex=0.6)
par(mfcol=c(2,2))
kjnj=log(da[,1])
djnj=diff(kjnj)
skjnj=diff(kjnj,4)
dkjnj=diff(skjnj)
acf(kjnj,lag=10)
acf(djnj,lag=10)
acf(skjnj,lag=10)
acf(dkjnj,lag=10)
c1=c("2","3","4","1")
c2=c("1","2","3","4")
par(mfcol=c(3,1))
plot(djnj,xlab='month',ylab='diff',type='l')
points(djnj,pch=c1,cex=0.7)
plot(skjnj,xlab='month',ylab='sea-diff',type='l')
points(skjnj,pch=c2,cex=0.7)
plot(dkjnj,xlab='month',ylab='dd',type='l')
points(dkjnj,pch=c1,cex=0.7)
# Estimation
m1=arima(kjnj,order=c(0,1,1),seasonal=list(order=c(0,1,1),period=4))
m1
##
## Call:
## arima(x = kjnj, order = c(0, 1, 1), seasonal = list(order = c(0, 1, 1), period = 4))
##
## Coefficients:
## ma1 sma1
## -1.0000 0.5436
## s.e. 0.0458 0.1147
##
## sigma^2 estimated as 0.3023: log likelihood = -62.35, aic = 130.7
tsdiag(m1,gof=20) # model checking
Box.test(m1$residuals,lag=12,type='Ljung')
##
## Box-Ljung test
##
## data: m1$residuals
## X-squared = 16.576, df = 12, p-value = 0.1662
pp=1-pchisq(13.30,10)
pp
## [1] 0.2073788
# Out-of-sample forecasting
kjnj=log(da[,3])
length(kjnj)
## [1] 78
y=kjnj[1:100]
m1=arima(y,order=c(0,1,1),seasonal=list(order=c(0,1,1),period=4))
m1
##
## Call:
## arima(x = y, order = c(0, 1, 1), seasonal = list(order = c(0, 1, 1), period = 4))
##
## Coefficients:
## ma1 sma1
## -0.0012 0.0035
## s.e. 0.0028 0.0021
##
## sigma^2 estimated as 1.476e-14: log likelihood = 1058.83, aic = -2111.67
# Prediction
pm1=predict(m1,10)
names(pm1)
## [1] "pred" "se"
pred=pm1$pred
##Exercise 8####################
da=read.table("q-GNPC96.txt",header=T)
head(da)
## year mon day gnp
## 1 1947 1 1 1780.4
## 2 1947 4 1 1778.1
## 3 1947 7 1 1776.6
## 4 1947 10 1 1804.0
## 5 1948 1 1 1833.4
## 6 1948 4 1 1867.6
gnp=log(da$gnp)
dgnp=diff(gnp)
m1=ar(dgnp,method='mle')
m1$order
## [1] 3
m2=arima(dgnp,order=c(3,0,0))
m2
##
## Call:
## arima(x = dgnp, order = c(3, 0, 0))
##
## Coefficients:
## ar1 ar2 ar3 intercept
## 0.3485 0.1386 -0.1317 0.0078
## s.e. 0.0616 0.0648 0.0617 0.0009
##
## sigma^2 estimated as 8.436e-05: log likelihood = 843.88, aic = -1677.76
m3=arima(dgnp,order=c(3,0,0),season=list(order=c(1,0,1),period=4))
m3
##
## Call:
## arima(x = dgnp, order = c(3, 0, 0), seasonal = list(order = c(1, 0, 1), period = 4))
##
## Coefficients:
## ar1 ar2 ar3 sar1 sma1 intercept
## 0.3385 0.1479 -0.1170 -0.5807 0.5294 0.0078
## s.e. 0.0627 0.0654 0.0638 0.4058 0.4203 0.0009
##
## sigma^2 estimated as 8.407e-05: log likelihood = 844.32, aic = -1674.64
source("backtest.R") # Perform backtest
mm2=backtest(m2,dgnp,215,1)
## [1] "RMSE of out-of-sample forecasts"
## [1] 0.007370678
## [1] "Mean absolute error of out-of-sample forecasts"
## [1] 0.004834576
mm3=backtest(m3,dgnp,215,1)
## [1] "RMSE of out-of-sample forecasts"
## [1] 0.007565077
## [1] "Mean absolute error of out-of-sample forecasts"
## [1] 0.005099862
mm1=ar(gnp,method='mle')
## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1
mm1$order # Find the identified order
## [1] 11
names(mm1)
## [1] "order" "ar" "var.pred" "x.mean"
## [5] "aic" "n.used" "n.obs" "order.max"
## [9] "partialacf" "resid" "method" "series"
## [13] "frequency" "call" "asy.var.coef"
print(mm1$aic,digits=3)
## 0 1 2 3 4 5 6 7 8
## 2129.29 137.23 89.26 6.80 6.17 7.52 12.66 5.92 5.69
## 9 10 11 12
## 5.10 6.41 0.00 2.60
aic=mm1$aic # For plotting below.
length(aic)
## [1] 13
plot(c(0:12),aic,type='h',xlab='order',ylab='aic')
lines(0:12,aic,lty=2)
da=read.table("q-GNPC96.txt",header=T)
head(da)
## year mon day gnp
## 1 1947 1 1 1780.4
## 2 1947 4 1 1778.1
## 3 1947 7 1 1776.6
## 4 1947 10 1 1804.0
## 5 1948 1 1 1833.4
## 6 1948 4 1 1867.6
G=da$gnp
LG=log(G)
LG
## [1] 7.484593 7.483301 7.482457 7.497762 7.513927 7.532409 7.537750
## [8] 7.539293 7.525101 7.521318 7.532195 7.522292 7.562058 7.592164
## [15] 7.631335 7.648644 7.660868 7.678234 7.698075 7.700159 7.709712
## [22] 7.710833 7.717351 7.749538 7.767857 7.775696 7.768998 7.753151
## [29] 7.748805 7.750055 7.761277 7.781723 7.810069 7.826244 7.839526
## [36] 7.845181 7.841375 7.849090 7.847997 7.863228 7.870433 7.868560
## [43] 7.877738 7.865687 7.838580 7.844789 7.867565 7.890395 7.910260
## [50] 7.935086 7.934334 7.938160 7.959975 7.955460 7.957352 7.944918
## [57] 7.951207 7.969219 7.985314 8.005467 8.023061 8.034631 8.043856
## [64] 8.047318 8.059814 8.071906 8.090556 8.098308 8.121005 8.132001
## [71] 8.145637 8.147838 8.172898 8.186520 8.205874 8.228871 8.253202
## [78] 8.256504 8.262869 8.271165 8.279799 8.279824 8.288384 8.295723
## [85] 8.316154 8.333006 8.340002 8.344244 8.359767 8.362409 8.368368
## [92] 8.363692 8.362292 8.364322 8.373023 8.361872 8.390019 8.396019
## [99] 8.403218 8.406329 8.424222 8.447393 8.457443 8.473617 8.499783
## [106] 8.511738 8.507890 8.517353 8.510430 8.512301 8.501450 8.496072
## [113] 8.482457 8.490069 8.507264 8.521703 8.543699 8.551479 8.556414
## [120] 8.563810 8.576612 8.595856 8.613285 8.611940 8.616749 8.653244
## [127] 8.663680 8.678053 8.679958 8.682165 8.691533 8.694134 8.697563
## [134] 8.676058 8.673445 8.689229 8.710257 8.701712 8.714288 8.703191
## [141] 8.686193 8.692926 8.687155 8.687543 8.699648 8.722303 8.741824
## [148] 8.762474 8.780665 8.797669 8.807158 8.814137 8.821732 8.830543
## [155] 8.844711 8.852879 8.861619 8.864040 8.874070 8.877452 8.883349
## [162] 8.894739 8.903217 8.920496 8.926690 8.938938 8.943584 8.957214
## [169] 8.966433 8.973681 8.981895 8.985107 8.995103 8.999199 8.998335
## [176] 8.993440 8.986534 8.991363 8.994991 9.000162 9.010950 9.021393
## [183] 9.031106 9.041661 9.044616 9.049937 9.056035 9.066816 9.077483
## [190] 9.090430 9.096500 9.107432 9.111194 9.113686 9.120120 9.128544
## [197] 9.136155 9.152086 9.160299 9.171745 9.178334 9.193947 9.205458
## [204] 9.212318 9.222733 9.231368 9.242953 9.260292 9.270400 9.278700
## [211] 9.290971 9.309787 9.311922 9.331504 9.331921 9.340210 9.335554
## [218] 9.342850 9.337880 9.347386 9.352126 9.355695 9.362091 9.364065
## [225] 9.366908 9.376541 9.392729 9.403932 9.411974 9.416142 9.424379
## [232] 9.430335 9.443308 9.446811 9.455394 9.458169 9.470757 9.474365
## [239] 9.473443 9.480925 9.482091 9.491632 9.503443 9.511259 9.505373
## [246] 9.508710 9.500574 9.472081 9.456075 9.454909 9.460702 9.470518
## [253] 9.479825 9.490877 9.497690 9.501808 9.505410 9.510756 9.515462
gnp=diff(LG)
dim(da)
## [1] 259 4
tdx=c(1:253)/4+1947 # create the time index
par(mfcol=c(2,1))
acf(gnp,lag=12)
pacf(gnp,lag=12) # compute PACF
m1=arima(gnp,order=c(3,0,0))
m1
##
## Call:
## arima(x = gnp, order = c(3, 0, 0))
##
## Coefficients:
## ar1 ar2 ar3 intercept
## 0.3485 0.1386 -0.1317 0.0078
## s.e. 0.0616 0.0648 0.0617 0.0009
##
## sigma^2 estimated as 8.436e-05: log likelihood = 843.88, aic = -1677.76