1.
theta_one<- theta$x
theta_two<- theta$y
theta_three<- theta$z
acf(theta_one)
acf(theta_two)
acf(theta_three)
pacf(theta_one)
pacf(theta_two)
pacf(theta_three)
b.
theta_three = ts(theta_three)
A = ts.intersect(theta_three, zL1 = lag(theta_three,-1), zL2 = lag(theta_three,-2), zL3 = lag(theta_three, -3), zL4 = lag(theta_three, -4), dframe = TRUE)
l1= lm(theta_three~ zL1 , data = A)
l2 = lm(theta_three~ zL2, data = A )
l3 = lm(theta_three~ zL3, data = A )
l4 = lm(theta_three~ zL4, data = A )
print(cbind(l1$coefficients, l2$coefficients, l3$coefficients, l4$coefficients))
## [,1] [,2] [,3] [,4]
## (Intercept) 0.006921002 0.01656603 0.02303718 0.02766529
## zL1 0.871979973 0.64303070 0.46110228 0.32084433
pacf4 = lm(theta_three~ zL1+zL2+zL3+zL4, data = A)
pacf4
##
## Call:
## lm(formula = theta_three ~ zL1 + zL2 + zL3 + zL4, data = A)
##
## Coefficients:
## (Intercept) zL1 zL2 zL3 zL4
## 0.006202 1.514386 -1.080288 0.637733 -0.239105
2.
a.
plot(rec)
b.
acf(rec)
pacf(rec)
c.
A1<- arima(rec, order = c(2,0,0))
A2<- arima(rec, order = c(1,0,0))
A1
##
## Call:
## arima(x = rec, order = c(2, 0, 0))
##
## Coefficients:
## ar1 ar2 intercept
## 1.3512 -0.4612 61.8585
## s.e. 0.0416 0.0417 4.0039
##
## sigma^2 estimated as 89.33: log likelihood = -1661.51, aic = 3331.02
A2
##
## Call:
## arima(x = rec, order = c(1, 0, 0))
##
## Coefficients:
## ar1 intercept
## 0.9249 61.2531
## s.e. 0.0177 6.4941
##
## sigma^2 estimated as 113.6: log likelihood = -1715.64, aic = 3437.27
d.
LinearTrendModel<- lm(rec~ time(rec))
summary(LinearTrendModel)
##
## Call:
## lm(formula = rec ~ time(rec))
##
## Residuals:
## Min 1Q Median 3Q Max
## -64.864 -22.696 5.365 23.036 45.855
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1035.5774 232.0641 -4.462 1.02e-05 ***
## time(rec) 0.5576 0.1179 4.731 3.00e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 27.34 on 451 degrees of freedom
## Multiple R-squared: 0.04728, Adjusted R-squared: 0.04517
## F-statistic: 22.38 on 1 and 451 DF, p-value: 2.996e-06
#significant at the five percent level
e.
ARTrend<- arima(rec, order = c(2,0,0), xreg = (time(rec)))
ARTrend
##
## Call:
## arima(x = rec, order = c(2, 0, 0), xreg = (time(rec)))
##
## Coefficients:
## ar1 ar2 intercept (time(rec))
## 1.3498 -0.4649 -846.4388 0.4614
## s.e. 0.0415 0.0417 680.8131 0.3458
##
## sigma^2 estimated as 89.01: log likelihood = -1660.67, aic = 3331.33
#not significant at the five percent level
f.
monthly<- ts(rec, frequency = 12)
Monthlyfactor <- factor(cycle(monthly))
MonthlyRegression<- lm(rec~ Monthlyfactor)
summary(MonthlyRegression)
##
## Call:
## lm(formula = rec ~ Monthlyfactor)
##
## Residuals:
## Min 1Q Median 3Q Max
## -61.515 -20.451 4.107 22.755 50.505
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 73.4708 4.2882 17.133 < 2e-16 ***
## Monthlyfactor2 -1.3182 6.0644 -0.217 0.828029
## Monthlyfactor3 -3.8024 6.0644 -0.627 0.530986
## Monthlyfactor4 -7.7679 6.0644 -1.281 0.200903
## Monthlyfactor5 -12.3763 6.0644 -2.041 0.041865 *
## Monthlyfactor6 -18.7755 6.0644 -3.096 0.002086 **
## Monthlyfactor7 -26.3671 6.0644 -4.348 1.71e-05 ***
## Monthlyfactor8 -29.2755 6.0644 -4.827 1.91e-06 ***
## Monthlyfactor9 -21.7911 6.0644 -3.593 0.000363 ***
## Monthlyfactor10 -8.6678 6.1053 -1.420 0.156392
## Monthlyfactor11 -4.0316 6.1053 -0.660 0.509374
## Monthlyfactor12 0.2341 6.1053 0.038 0.969434
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 26.43 on 441 degrees of freedom
## Multiple R-squared: 0.129, Adjusted R-squared: 0.1073
## F-statistic: 5.937 on 11 and 441 DF, p-value: 4.848e-09
AIC(MonthlyRegression)
## [1] 4266.236
#this model does not do as well, Higher AIC than other models.
g.
ARMA<- dynlm(rec~ lag(rec, -1)+ lag(rec, -2)+ Monthlyfactor)
AIC(ARMA)
## [1] 3248.416
plot(ARMA$residuals)
acf(ARMA$residuals, lag.max = 50)
pacf(ARMA$residuals, lag.max = 50)
3.
a.
rec3<- window(rec, end=c(1985, 8))
ActualRec<- window(rec, start=c(1985,8), end=c(1987,7))
ARTrend3<- arima(rec3, order = c(2,0,0), include.mean = FALSE)
prediction<- forecast(ARTrend3, h=24)[4]
plot(forecast(ARTrend3, h=24))
predict<- unlist(prediction)
RMSE <- sqrt(mean((ActualRec-predict)^2))
RMSE
## [1] 25.90144
b.
rec4 <- rec[1:428]
dd<- ts(rec4, frequency = 12)
Monthlyfactor4 <- factor(cycle(dd))
Dummy1<- model.matrix(rec4~ Monthlyfactor4)
Dummy1<- Dummy1[,-1]
Monthlyfactor44 <- factor(cycle(ActualRec))
Dummy4<- model.matrix(ActualRec~ Monthlyfactor44)
Dummy4<- Dummy4[,-1]
ARTrend3b<- arima(rec4, order = c(2,0,0), include.mean = FALSE, xreg = Dummy1)
AIC(ARTrend3b)
## [1] 3112.305
plot(ARTrend3b)
predict3b<- predict(ARTrend3b, n.ahead=24, newxreg = Dummy4)
rec4<- ts(rec4)
ts.plot(rec4, predict3b$pred, col=c(1,2))
actual<- ts(ActualRec, frequency = 12)
pred<- ts(predict3b$pred, frequency = 12)
MEAN<- mean(actual- pred)
RMSE2<- rmse(actual, pred)
RMSE2
## [1] 12.11171