library(quantmod)
## Loading required package: xts
## Loading required package: zoo
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
## Loading required package: TTR
## Version 0.4-0 included new data defaults. See ?getSymbols.
library(fUnitRoots)
## Loading required package: timeDate
## Loading required package: timeSeries
##
## Attaching package: 'timeSeries'
## The following object is masked from 'package:zoo':
##
## time<-
## Loading required package: fBasics
##
## Attaching package: 'fBasics'
## The following object is masked from 'package:TTR':
##
## volatility
library(timeSeries)
library(fBasics)
data1=read.table("m-unrate-4811.txt", header=T)
#=============================
#Exersice1
#=============================
head(data)
##
## 1 function (..., list = character(), package = NULL, lib.loc = NULL,
## 2 verbose = getOption("verbose"), envir = .GlobalEnv)
## 3 {
## 4 fileExt <- function(x) {
## 5 db <- grepl("\\\\.[^.]+\\\\.(gz|bz2|xz)$", x)
## 6 ans <- sub(".*\\\\.", "", x)
values = seq(from = as.Date("1948-01-01"), to = as.Date("2011-11-03"), by = 'month')
rownames(data1) <- values
str(data1)
## 'data.frame': 767 obs. of 1 variable:
## $ rate: num 3.4 3.8 4 3.9 3.5 3.6 3.6 3.9 3.8 3.7 ...
gdp=log(data1[,1])
m1=ar(diff(gdp),method="mle")
m1$order
## [1] 12
adfTest(gdp,lags=10,type=c("c"))
##
## Title:
## Augmented Dickey-Fuller Test
##
## Test Results:
## PARAMETER:
## Lag Order: 10
## STATISTIC:
## Dickey-Fuller: -3.3319
## P VALUE:
## 0.01523
##
## Description:
## Thu Oct 17 23:00:44 2019 by user: Minjin
eps=log(data1$rate)
koeps=ts(eps,frequency=4,start=c(1948,1))
plot(koeps,type="l")
points(koeps,pch=0,cex=0.6)
#Obtain ACF plot
par(mfcol=c(2,2))
koeps=log(data1$rate)
deps=diff(koeps)
sdeps=diff(koeps,4)
ddeps=diff(sdeps)
acf(koeps,lag=20)
acf(deps,lag=20)
acf(sdeps,lag=20)
acf(ddeps,lag=20)
#Obtain time plots
c1=c("2","3","4","1")
c2=c("1","2","3","4")
par(mfcol=c(3,1))
plot(deps,xlab="year",ylab="diff",type="l")
points(deps,pch=c1,cex=0.7)
plot(sdeps,xlab="year",ylab="sea-diff",type="l")
points(sdeps,pch=c2,cex=0.7)
plot(ddeps,xlab="year",ylab="dd",type="l")
points(ddeps,pch=c1,cex=0.7)
#Estimation
m1=arima(koeps,order=c(0,1,1),seasonal=list(order=c(0,1,1),period=4))
m1
##
## Call:
## arima(x = koeps, 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
koeps=log(data1$rate)
length(koeps)
## [1] 767
y=koeps[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,7)
names(pm1)
## [1] "pred" "se"
pred=pm1$pred
se=pm1$se
#exercise 7
#===================================
data2=read.table("q-jnj-earns-9211.txt", header=T)
head(data2)
## 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
J=data2$earns
JN=ts(J,frequency=12,start=c(1992,1))
par(mfcol=c(2,1))
plot(JN,xlab="year",ylab="returns")
title(main="(a): Simple returns")
acf(J,lag=24)
ln.J=log(J+1)
Box.test(J,lag=12,type="Ljung")
##
## Box-Ljung test
##
## data: J
## X-squared = 584.59, df = 12, p-value < 2.2e-16
Box.test(ln.J,lag=12,type="Ljung")
##
## Box-Ljung test
##
## data: ln.J
## X-squared = 598.48, df = 12, p-value < 2.2e-16
JNJ=ts(J,frequency=12,start=c(1992,1))
par(mfcol=c(2,1))
plot(J,xlab="year",ylab="returns")
title(main="(a): Simple returns")
acf(JN,lag=24)
gnp=diff(ln.J)
dim(data2)
## [1] 78 4
tdx=c(1:78)/4+1992
par(mfcol=c(2,1))
plot(tdx,ln.J,xlab="year",ylab="gnp",type="l")
plot(tdx[2:78],gnp,type="l",xlab="year",ylab="growth")
acf(gnp,lag=12)
pacf(gnp,lag=12)
m1=arima(gnp,order=c(3,0,0))
m1
##
## Call:
## arima(x = gnp, order = c(3, 0, 0))
##
## Coefficients:
## ar1 ar2 ar3 intercept
## -0.9547 -0.9342 -0.9461 0.0087
## s.e. 0.0324 0.0363 0.0277 0.0006
##
## sigma^2 estimated as 0.0003623: log likelihood = 192.15, aic = -374.3
tsdiag(m1,gof=12)
p1=c(1,-m1$coef[1:3])
r1=polyroot(p1)
r1
## [1] 0.014824+1.019333i -1.017089+0.000000i 0.014824-1.019333i
Mod(r1)
## [1] 1.019441 1.017089 1.019441
k=2*pi/acos(0.014824/1.019441)
k
## [1] 4.037376
mm1=ar(gnp,method="mle")
mm1$order
## [1] 4
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
## 251.685 232.820 226.296 51.073 0.000 1.821 3.812 5.118 0.102
## 9 10 11 12
## 1.430 1.319 2.384 4.215
aic=mm1$aic
length(aic)
## [1] 13
plot(c(0:12),aic,type="h",xlab="order",ylab="aic")
lines(0:12,aic,lty=2)
m1=arima(J,order=c(0,0,9))
m1
##
## Call:
## arima(x = J, order = c(0, 0, 9))
##
## Coefficients:
## ma1 ma2 ma3 ma4 ma5 ma6 ma7 ma8
## 1.0058 1.2457 1.2014 2.1391 1.5133 1.319 1.0087 1.2378
## s.e. 0.1272 0.1732 0.2721 0.2876 0.2363 0.215 0.2427 0.2446
## ma9 intercept
## 0.3804 0.6094
## s.e. 0.1330 0.0791
##
## sigma^2 estimated as 0.003708: log likelihood = 96.24, aic = -170.47
m1=arima(J,order=c(0,0,9),fixed=c(NA,0,NA,0,0,0,0,0,NA,NA))
m1
##
## Call:
## arima(x = J, 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
## 0.1294 0 0.6503 0 0 0 0 0 -0.2064 0.6002
## s.e. 0.2891 0 0.3004 0 0 0 0 0 0.2011 0.0468
##
## sigma^2 estimated as 0.06913: log likelihood = -7.43, aic = 24.87
sqrt(0.06913)
## [1] 0.2629258
Box.test(m1$residuals,lag=12,type="Ljung")
##
## Box-Ljung test
##
## data: m1$residuals
## X-squared = 394.68, df = 12, p-value < 2.2e-16
pv=1-pchisq(394.68,9)
pv
## [1] 0
pv=1-pchisq(394.68,9)
pv
## [1] 0
m1=arima(J[1:78],order=c(0,0,9),fixed=c(NA,0,NA,0,0,0,0,0,NA,NA))
m1
##
## Call:
## arima(x = J[1:78], 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
## 0.1294 0 0.6503 0 0 0 0 0 -0.2064 0.6002
## s.e. 0.2891 0 0.3004 0 0 0 0 0 0.2011 0.0468
##
## sigma^2 estimated as 0.06913: log likelihood = -7.43, aic = 24.87
predict(m1,10)
## $pred
## Time Series:
## Start = 79
## End = 88
## Frequency = 1
## [1] 0.9486664 0.6352115 0.8469624 0.5757401 0.4878375 0.5478441 0.4812723
## [8] 0.5757913 0.4822479 0.6002087
##
## $se
## Time Series:
## Start = 79
## End = 88
## Frequency = 1
## [1] 0.2629332 0.2651269 0.2651269 0.3154777 0.3154777 0.3154777 0.3154777
## [8] 0.3154777 0.3154777 0.3201099
# exercise 8
#========================
library(fBasics)
data3=read.table("q-GNPC96.txt",header=T)
head(data3)
## 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
gdp=log(data3$gnp)
dgdp=diff(gdp)
m1=ar(dgdp,method="mle")
m1$order
## [1] 3
m2=arima(dgdp,order=c(4,0,0))
m2
##
## Call:
## arima(x = dgdp, order = c(4, 0, 0))
##
## Coefficients:
## ar1 ar2 ar3 ar4 intercept
## 0.3369 0.1513 -0.1010 -0.0887 0.0078
## s.e. 0.0619 0.0652 0.0651 0.0619 0.0008
##
## sigma^2 estimated as 8.368e-05: log likelihood = 844.9, aic = -1677.8
m3=arima(dgdp,order=c(3,0,0),season=list(order=c(1,0,1),period=4))
m3
##
## Call:
## arima(x = dgdp, 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,dgdp,215,1)
## [1] "RMSE of out-of-sample forecasts"
## [1] 0.007515004
## [1] "Mean absolute error of out-of-sample forecasts"
## [1] 0.005053905
mm3=backtest(m3,dgdp,215,1)
## [1] "RMSE of out-of-sample forecasts"
## [1] 0.007565077
## [1] "Mean absolute error of out-of-sample forecasts"
## [1] 0.005099862