Task 8 Code and Output

library(dynlm)

## Loading required package: zoo

##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## [Link], [Link]
library(ggplot2)
library(AER)

## Loading required package: car

## Loading required package: carData
## Loading required package: lmtest
## Loading required package: sandwich
## Loading required package: survival
library(Hmisc)

## Loading required package: lattice

## Loading required package: Formula
##
## Attaching package: 'Hmisc'
## The following objects are masked from 'package:base':
##
## [Link], units
library(forecast)

## Registered S3 method overwritten by 'quantmod':

## method from
## [Link] zoo
library(x12)

## Loading required package: x13binary

## x12 is ready to use.
## Use the package x12GUI for a Graphical User Interface.
## By default the X13-ARIMA-SEATS binaries provided by the R package x13binary
## are used but this can be changed with x12path(validpath)
## ---------------

1
## Suggestions and bug-reports can be submitted at: [Link]
library(dLagM)

## Loading required package: nardl

##
## Attaching package: 'dLagM'
## The following object is masked from 'package:forecast':
##
## forecast
library(TSA)

## Registered S3 methods overwritten by 'TSA':

## method from
## [Link] forecast
## [Link] forecast
##
## Attaching package: 'TSA'
## The following objects are masked from 'package:stats':
##
## acf, arima
## The following object is masked from 'package:utils':
##
## tar
# TASK 1
aerosol = [Link]("[Link]")

class(aerosol)

## [1] "[Link]"
# Be careful about the transpose t() operation here!
aerosol = ts([Link](t([Link](aerosol[,2:13]))),
start=c(1986,1),frequency = 12)
class(aerosol)

## [1] "ts"
head(aerosol)

## Jan Feb Mar Apr May Jun

## 1986 0.07644696 0.06693375 0.05330857 0.03857360 0.03472612 0.01393000
par(mfrow=c(1,1))
plot(aerosol, ylab="Monthly median aerosol optical depth",xlab="Year",
main="Time series plot of mean monthly median aerosol
optical depth at Cape Grim")
points(y=aerosol,x=time(aerosol), pch=[Link](season(aerosol)))

2
 Time series plot of mean monthly median aerosol
optical depth at Cape Grim
Monthly median aerosol optical depth

0.20

S
JN
S O
N
DAJ
0.15

O
MD
AJFA
M FJ F
0.10

JJANJ J
AS
M D N AD N
O D O
MA AM D S D JF N
M J FJ N F J
J S
F SF MO DA
M NA
JSM OJM
JA
S F JMMD D
FA JM NS A MJA
SMJS
O
NJ D
MA
A
F
O
SN JA MONM F S M S O
SM MO A
J J
N
D
OJ O N OJ A J S
0.05

OF DJ O D S FSJ MA JA S SD
M O MNM JA
S
MM JJNJJ J OM OMJSMMN
A
FS AN
M
D
NJ
M A
DF
MJ
OANJA J OAS FJ
A
DJDM A NF S
J
O
JAA F N A M N
A D N DFA J
M
F J JMOAJJ
MJA
A J
J J JA M NJAJ DF M
A OS F FA
J
O
D
S
O JAN JM M FS MNM
S
OA J
MO DJD J NMD
M N FA
O
J MJADFMA
FM J M JA J DM
F MAA
AA A A J J MA J JAD M
J J J J M JM J JA A J

1985 1990 1995 2000 2005 2010 2015

Year

# par(mfrow=c(1,2)) # Put the ACF and PACF plots next to each other
acf(aerosol, [Link] = 48, main = "Sample ACF for the median monthly aerosol series")

3
 Sample ACF for the median monthly aerosol series
0.8
0.6
0.4
ACF

0.2
0.0

0 1 2 3 4

Lag

# pacf(aerosol, [Link] = 48, main = "Sample PACF for the median monthly aerosol series")

aerosol.x12 = x12(aerosol)
plot(aerosol.x12 , sa=TRUE , trend=TRUE)

4
 0.20
0.15
Original Series, Seasonally Adjusted Series and Trend
Value

0.10
0.05

1986.1 1990.1 1995.1 2000.1 2005.1 2010.1 2015.1

Date

Original Seasonally Adjusted Trend

fit.AAA1 = ets(aerosol, model = "AAA", damped = FALSE, bounds="admiss" )

summary(fit.AAA1)

## ETS(A,A,A)
##
## Call:
## ets(y = aerosol, model = "AAA", damped = FALSE, bounds = "admiss")
##
## Smoothing parameters:
## alpha = 0.077
## beta = 0
## gamma = 0.1077
##
## Initial states:
## l = 0.0087
## b = 0
## s = -0.0108 0.0171 0.0191 0.0272 -0.0152 -0.0253
## -0.0075 -0.0408 -0.0143 -0.0075 0.0276 0.0303
##
## sigma: 0.0247
##
## AIC AICc BIC
## -520.8898 -519.0352 -455.4023
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE

5
## Training set 0.0007434413 0.02415276 0.01780343 -12.85203 39.56627 0.7933962
## ACF1
## Training set 0.5480029
checkresiduals(fit.AAA1)

Residuals from ETS(A,A,A)

0.10

0.05

0.00

−0.05

1985 1990 1995 2000 2005 2010 2015

50

0.4 40
count

30
ACF

0.2
20

0.0 10

0
12 24 36 −0.05 0.00 0.05 0.10
Lag residuals

##
## Ljung-Box test
##
## data: Residuals from ETS(A,A,A)
## Q* = 345.1, df = 8, p-value < 2.2e-16
##
## Model df: 16. Total lags used: 24
fit.AAA2 = ets(aerosol, damped = FALSE, model = "AAA")
summary(fit.AAA2)

## ETS(A,A,A)
##
## Call:
## ets(y = aerosol, model = "AAA", damped = FALSE)
##
## Smoothing parameters:
## alpha = 0.2814
## beta = 7e-04
## gamma = 0.0436
##
## Initial states:

6
## l = -0.0238
## b = -1e-04
## s = 0.0082 0.0074 0.0105 0.0147 -0.0023 -0.0136
## -0.0168 -0.0123 -0.0041 -6e-04 -0.007 0.0161
##
## sigma: 0.0201
##
## AIC AICc BIC
## -665.0989 -663.2444 -599.6115
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.0004839745 0.0196328 0.01416031 -10.52223 31.4338 0.6310433
## ACF1
## Training set 0.2976457
checkresiduals(fit.AAA2)

Residuals from ETS(A,A,A)

0.08

0.04

0.00

−0.04

1985 1990 1995 2000 2005 2010 2015

0.3

0.2
40
count
ACF

0.1

20
0.0

−0.1 0
12 24 36 −0.05 0.00 0.05
Lag residuals

##
## Ljung-Box test
##
## data: Residuals from ETS(A,A,A)
## Q* = 62.587, df = 8, p-value = 1.445e-10
##
## Model df: 16. Total lags used: 24

7
fit.AAA3 = ets(aerosol, damped = FALSE, model = "AAA", upper=rep(1,4))
summary(fit.AAA3)

## ETS(A,A,A)
##
## Call:
## ets(y = aerosol, model = "AAA", damped = FALSE, upper = rep(1,
##
## Call:
## 4))
##
## Smoothing parameters:
## alpha = 0.2814
## beta = 7e-04
## gamma = 0.0436
##
## Initial states:
## l = -0.0238
## b = -1e-04
## s = 0.0082 0.0074 0.0105 0.0147 -0.0023 -0.0136
## -0.0168 -0.0123 -0.0041 -6e-04 -0.007 0.0161
##
## sigma: 0.0201
##
## AIC AICc BIC
## -665.0992 -663.2446 -599.6117
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.0004838914 0.01963279 0.01416031 -10.52239 31.43382 0.6310432
## ACF1
## Training set 0.2976446
checkresiduals(fit.AAA3)

8
 Residuals from ETS(A,A,A)
0.08

0.04

0.00

−0.04

1985 1990 1995 2000 2005 2010 2015

0.3

0.2
40

count
ACF

0.1

20
0.0

−0.1 0
12 24 36 −0.05 0.00 0.05
Lag residuals

##
## Ljung-Box test
##
## data: Residuals from ETS(A,A,A)
## Q* = 62.587, df = 8, p-value = 1.445e-10
##
## Model df: 16. Total lags used: 24
fit.AAA4 = ets(aerosol, damped = FALSE, model = "AAA", [Link] = "mse")
summary(fit.AAA4)

## ETS(A,A,A)
##
## Call:
## ets(y = aerosol, model = "AAA", damped = FALSE, [Link] = "mse")
##
## Smoothing parameters:
## alpha = 0.2814
## beta = 7e-04
## gamma = 0.0436
##
## Initial states:
## l = -0.0238
## b = -1e-04
## s = 0.0082 0.0074 0.0105 0.0147 -0.0023 -0.0136
## -0.0168 -0.0123 -0.0041 -6e-04 -0.007 0.0161
##

9
## sigma: 0.0201
##
## AIC AICc BIC
## -665.0989 -663.2444 -599.6115
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.0004839745 0.0196328 0.01416031 -10.52223 31.4338 0.6310433
## ACF1
## Training set 0.2976457
checkresiduals(fit.AAA4)

Residuals from ETS(A,A,A)

0.08

0.04

0.00

−0.04

1985 1990 1995 2000 2005 2010 2015

0.3

0.2
40
count
ACF

0.1

20
0.0

−0.1 0
12 24 36 −0.05 0.00 0.05
Lag residuals

##
## Ljung-Box test
##
## data: Residuals from ETS(A,A,A)
## Q* = 62.587, df = 8, p-value = 1.445e-10
##
## Model df: 16. Total lags used: 24
fit.AAA5 = ets(aerosol, damped = FALSE, model = "AAA", [Link] = "mae")
summary(fit.AAA5)

## ETS(A,A,A)
##
## Call:
## ets(y = aerosol, model = "AAA", damped = FALSE, [Link] = "mae")

10
##
## Smoothing parameters:
## alpha = 0.0423
## beta = 1e-04
## gamma = 0.0379
##
## Initial states:
## l = 0.0366
## b = -3e-04
## s = 0.0086 0.0031 -0.0072 -0.0048 0.0051 -0.0108
## -0.0068 -0.0046 5e-04 0.0074 0.0038 0.0057
##
## sigma: 0.0268
##
## AIC AICc BIC
## -465.1170 -463.2625 -399.6296
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.004564041 0.02616786 0.01734157 -9.107762 36.10969 0.7728136
## ACF1
## Training set 0.6755666
checkresiduals(fit.AAA5)

Residuals from ETS(A,A,A)

0.10

0.05

0.00

−0.05
1985 1990 1995 2000 2005 2010 2015

0.50 40
count
ACF

0.25
20

0.00
0
12 24 36 −0.05 0.00 0.05 0.10
Lag residuals

##
## Ljung-Box test

11
##
## data: Residuals from ETS(A,A,A)
## Q* = 781.45, df = 8, p-value < 2.2e-16
##
## Model df: 16. Total lags used: 24
fit.AAdA1 = ets(aerosol, model = "AAA", damped = TRUE, bounds="admiss" )
summary(fit.AAdA1)

## ETS(A,Ad,A)
##
## Call:
## ets(y = aerosol, model = "AAA", damped = TRUE, bounds = "admiss")
##
## Smoothing parameters:
## alpha = 0.3774
## beta = -0.0171
## gamma = 0.0455
## phi = 0.9496
##
## Initial states:
## l = 0.0629
## b = 5e-04
## s = 0.0064 0.0157 0.0104 0.0139 -0.0038 -0.017
## -0.0067 -0.0144 -0.0103 -0.0011 -0.0018 0.0088
##
## sigma: 0.0185
##
## AIC AICc BIC
## -722.0682 -719.9892 -652.7285
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -0.001063393 0.01803791 0.01342327 -15.4257 31.89714 0.5981976
## ACF1
## Training set 0.166058
checkresiduals(fit.AAdA1)

12
 Residuals from ETS(A,Ad,A)

0.05

0.00

−0.05
1985 1990 1995 2000 2005 2010 2015

60

0.1
40

count
ACF

0.0
20

−0.1
0
12 24 36 −0.05 0.00 0.05
Lag residuals

##
## Ljung-Box test
##
## data: Residuals from ETS(A,Ad,A)
## Q* = 42.51, df = 7, p-value = 4.147e-07
##
## Model df: 17. Total lags used: 24
fit.AAdA2 = ets(aerosol, model = "AAA", damped = TRUE)
summary(fit.AAdA2)

## ETS(A,Ad,A)
##
## Call:
## ets(y = aerosol, model = "AAA", damped = TRUE)
##
## Smoothing parameters:
## alpha = 0.2307
## beta = 2e-04
## gamma = 1e-04
## phi = 0.9769
##
## Initial states:
## l = 0.0706
## b = -7e-04
## s = 0.0043 0.0075 0.0062 0.0165 -0.0022 -0.0114
## -0.0083 -0.008 -0.0088 -0.0055 -1e-04 0.0098

13
##
## sigma: 0.0193
##
## AIC AICc BIC
## -693.9111 -691.8321 -624.5715
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 6.016721e-05 0.01878261 0.01377787 -12.89539 31.63064 0.6140001
## ACF1
## Training set 0.3466
checkresiduals(fit.AAdA2)

Residuals from ETS(A,Ad,A)

0.05

0.00

−0.05

1985 1990 1995 2000 2005 2010 2015

0.3

0.2 40
count
ACF

0.1
20
0.0

−0.1
0
12 24 36 −0.05 0.00 0.05
Lag residuals

##
## Ljung-Box test
##
## data: Residuals from ETS(A,Ad,A)
## Q* = 88.559, df = 7, p-value = 2.22e-16
##
## Model df: 17. Total lags used: 24
fit.AAdA3 = ets(aerosol, model = "AAA", damped = TRUE, upper=rep(1,4))
summary(fit.AAdA3)

## ETS(A,Ad,A)
##
## Call:

14
## ets(y = aerosol, model = "AAA", damped = TRUE, upper = rep(1,
##
## Call:
## 4))
##
## Smoothing parameters:
## alpha = 0.3027
## beta = 0.0297
## gamma = 0.0062
## phi = 0.9979
##
## Initial states:
## l = 0.0397
## b = -0.0018
## s = 0.0014 0.014 0.0083 0.0124 -0.0042 -0.0165
## -0.013 -0.0101 -0.0037 -6e-04 0.0013 0.0108
##
## sigma: 0.0193
##
## AIC AICc BIC
## -691.5071 -689.4280 -622.1674
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.0002122696 0.01884759 0.01393632 -9.117744 29.93654 0.6210615
## ACF1
## Training set 0.2647326
checkresiduals(fit.AAdA3)

15
 Residuals from ETS(A,Ad,A)

0.05

0.00

−0.05

1985 1990 1995 2000 2005 2010 2015

60
0.2

0.1 40

count
ACF

0.0
20

−0.1
0
12 24 36 −0.05 0.00 0.05
Lag residuals

##
## Ljung-Box test
##
## data: Residuals from ETS(A,Ad,A)
## Q* = 67.266, df = 7, p-value = 5.263e-12
##
## Model df: 17. Total lags used: 24
fit.AAdA4 = ets(aerosol, model = "AAA", damped = TRUE, [Link] = "mse")
summary(fit.AAdA4)

## ETS(A,Ad,A)
##
## Call:
## ets(y = aerosol, model = "AAA", damped = TRUE, [Link] = "mse")
##
## Smoothing parameters:
## alpha = 0.2307
## beta = 2e-04
## gamma = 1e-04
## phi = 0.9769
##
## Initial states:
## l = 0.0706
## b = -7e-04
## s = 0.0043 0.0075 0.0062 0.0165 -0.0022 -0.0114
## -0.0083 -0.008 -0.0088 -0.0055 -1e-04 0.0098

16
##
## sigma: 0.0193
##
## AIC AICc BIC
## -693.9111 -691.8321 -624.5715
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 6.016721e-05 0.01878261 0.01377787 -12.89539 31.63064 0.6140001
## ACF1
## Training set 0.3466
checkresiduals(fit.AAdA4)

Residuals from ETS(A,Ad,A)

0.05

0.00

−0.05

1985 1990 1995 2000 2005 2010 2015

0.3

0.2 40
count
ACF

0.1
20
0.0

−0.1
0
12 24 36 −0.05 0.00 0.05
Lag residuals

##
## Ljung-Box test
##
## data: Residuals from ETS(A,Ad,A)
## Q* = 88.559, df = 7, p-value = 2.22e-16
##
## Model df: 17. Total lags used: 24
fit.AAdA5 = ets(aerosol, model = "AAA", damped = TRUE, [Link] = "mae")
summary(fit.AAdA5)

## ETS(A,Ad,A)
##
## Call:

17
## ets(y = aerosol, model = "AAA", damped = TRUE, [Link] = "mae")
##
## Smoothing parameters:
## alpha = 0.2242
## beta = 0.0305
## gamma = 1e-04
## phi = 0.9606
##
## Initial states:
## l = 0.0645
## b = -0.0031
## s = 0.0047 0.0062 0.0048 0.0094 -0.0037 -0.0133
## -0.0094 -0.0076 -0.0043 0.0016 0.0037 0.0079
##
## sigma: 0.0194
##
## AIC AICc BIC
## -688.1761 -686.0971 -618.8365
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.0002455318 0.01893801 0.0135568 -10.04702 29.61988 0.6041483
## ACF1
## Training set 0.3421299
checkresiduals(fit.AAdA5)

Residuals from ETS(A,Ad,A)

0.05

0.00

−0.05

1985 1990 1995 2000 2005 2010 2015

60
0.3

0.2 40
count
ACF

0.1

0.0 20

−0.1
0
12 24 36 −0.05 0.00 0.05
Lag residuals

18
##
## Ljung-Box test
##
## data: Residuals from ETS(A,Ad,A)
## Q* = 96.175, df = 7, p-value < 2.2e-16
##
## Model df: 17. Total lags used: 24
# Task 2
A = ts(matrix(NA,120,5000),start=c(2015,1),frequency = 12)
# n = length(aerosol)
# h = 10
M = 5000
for (i in 1:M){
A[,i] = simulate(fit.AAdA1 , initstate = fit.AAdA1 $states[25,] , nsim=120)
# Generate random epsilons and apply model formulation
}

par(mfrow=c(1,1))
plot(aerosol , ylim=range(aerosol,A) , xlim=c(1985,2026) ,
ylab="Median aerosol optical depth" , xlab="Year")
for(i in 1:10){
lines(A[,i],col="gray")
}
text(1990,0,"Historical data",adj=0)
text(2015,0,"10 simulated future sample paths",adj=0)
0.20
Median aerosol optical depth

0.10
0.00

Historical data 10 simulated future sample paths

−0.10

1990 2000 2010 2020

Year

19
plot(aerosol , ylim=range(aerosol,A) , xlim=c(1985,2026) ,
ylab="Median aerosol optical depth" , xlab="Year")
for(i in 1:20){
lines(A[,i],col="gray")
}
text(1990,0,"Historical data",adj=0)
text(2015,0,"20 simulated future sample paths",adj=0)
0.20
Median aerosol optical depth

0.10
0.00

Historical data 20 simulated future sample paths

−0.10

1990 2000 2010 2020

Year

plot(aerosol , ylim=range(aerosol,A) , xlim=c(1985,2026) ,

ylab="Median aerosol optical depth" , xlab="Year")
for(i in 1:30){
lines(A[,i],col="gray")
}
text(1990,0,"Historical data",adj=0)
text(2015,0,"30 simulated future sample paths",adj=0)

20
 0.20
Median aerosol optical depth

0.10
0.00

Historical data 30 simulated future sample paths

−0.10

1990 2000 2010 2020

Year

plot(aerosol , ylim=range(aerosol,A) , xlim=c(1985,2026) ,

ylab="Median aerosol optical depth" , xlab="Year")
for(i in 1:100){
lines(A[,i],col="gray")
}
text(1990,0,"Historical data",adj=0)
text(2015,0,"100 simulated future sample paths",adj=0)

21
 0.20
Median aerosol optical depth

0.10
0.00

Historical data 100 simulated future sample paths

−0.10

1990 2000 2010 2020

Year

plot(aerosol , ylim=range(aerosol,A) , xlim=c(1985,2026) ,

ylab="Median aerosol optical depth" , xlab="Year")
for(i in 1:1000){
lines(A[,i],col="gray")
}
text(1990,0,"Historical data",adj=0)
text(2015,0,"5000 simulated future sample paths",adj=0)

22
 0.20
Median aerosol optical depth

0.10
0.00

Historical data 5000 simulated future sample paths

−0.10

1990 2000 2010 2020

Year

N = 120
data = aerosol
xlim=c(1985,2026)
Pi = array(NA, dim=c(N,2))
avrg = array(NA, N)
# Calcualte the interval estimates and mid point
for (i in 1:N){
Pi[i,] = quantile(A[i,],type=8,prob=c(.05,.95))
avrg[i] = mean(A[i,]) # This would be median as well
}
# Create ts objects for plotting
[Link] = ts(Pi[,1],start=end(data),f=12)
[Link] = ts(Pi[,2],start=end(data),f=12)
[Link] = ts(avrg,start=end(data),f=12)

plot(data,xlim=xlim , ylim=range(data,A),ylab="Y",xlab="Year", main="

10 years ahead predictions")
lines([Link],col="blue", type="l")
lines([Link],col="red", type="l")
lines([Link],col="green", type="l")
legend("topleft", lty=1, pch=1, col=c("black","blue","red","green"), [Link] = 4,
c("Data","5% lower limit","95% upper limit","Mean prediction"))

23
 0.20 10 years ahead predictions

Data
5% lower limit
95% upper limit
Mean prediction
0.10
Y

0.00
−0.10

1990 2000 2010 2020

Year

24