This is using Model 1 (direct fit with ARMA) from the last lecture.
library(forecast)
source('https://nmimoto.github.io/R/TS-00.txt')
= read.csv("https://nmimoto.github.io/datasets/lake.csv")
D = ts(D[,1], start=1875, freq=1)
Lake length(Lake)
## [1] 97
plot(Lake, type="o")
= auto.arima(Lake, d=0, stepwise=FALSE, approximation=FALSE)
Fit1 Fit1
## Series: Lake
## ARIMA(1,0,1) with non-zero mean
##
## Coefficients:
## ar1 ma1 mean
## 0.7665 0.3393 9.1290
## s.e. 0.0773 0.1123 0.3861
##
## sigma^2 = 0.4784: log likelihood = -101.09
## AIC=210.18 AICc=210.62 BIC=220.48
Randomness.tests(Fit1$residuals)
## B-L test H0: the series is uncorrelated
## M-L test H0: the square of the series is uncorrelated
## J-B test H0: the series came from Normal distribution
## SD : Standard Deviation of the series
## BL15 BL20 BL25 ML15 ML20 JB SD
## [1,] 0.963 0.952 0.934 0.567 0.641 0.89 0.684
plot(forecast(Fit1, 10))
abline(h=mean(Lake), col="red")
forecast(Fit1, 10)
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 1972 9.752287 8.865860 10.63871 8.396614 11.10796
## 1973 9.606775 8.285205 10.92834 7.585609 11.62794
## 1974 9.495234 7.975010 11.01546 7.170252 11.82022
## 1975 9.409735 7.784069 11.03540 6.923494 11.89598
## 1976 9.344197 7.659652 11.02874 6.767908 11.92049
## 1977 9.293960 7.575760 11.01216 6.666200 11.92172
## 1978 9.255451 7.517780 10.99312 6.597913 11.91299
## 1979 9.225933 7.476923 10.97494 6.551053 11.90081
## 1980 9.203307 7.447668 10.95895 6.518289 11.88832
## 1981 9.185963 7.426440 10.94548 6.495006 11.87692
source('https://nmimoto.github.io/R/TS-00.txt')
= Rolling1step.forecast(Lake,
Pred window.size=60,
Arima.order=c(1,0,1),
include.mean=TRUE)
##
## Total length 97 , window size 60 .
## Last 37 obs retrospectively forecasted with Rolling 1-step
## prediction using same order and fized window size.
##
## Average Prediction Error: 0.0058
## root Mean Squared Error of Prediction: 0.7303
# Pred
This is using Model 2 (linear trend + ARMA).
= Arima(Lake, order=c(1,0,1), xreg=time(Lake) )
Fit2 Fit2
## Series: Lake
## Regression with ARIMA(1,0,1) errors
##
## Coefficients:
## ar1 ma1 intercept xreg
## 0.6682 0.3817 55.5443 -0.0242
## s.e. 0.0936 0.1136 18.0324 0.0094
##
## sigma^2 = 0.4612: log likelihood = -98.66
## AIC=207.33 AICc=207.98 BIC=220.2
=lm(Lake~time(Lake))
Reg2=10
hplot(forecast(Fit2, h, xreg=last(time(Lake))+(1:h)/frequency(Lake)))
abline(Reg2, col="red")
forecast(Fit2, h, xreg=last(time(Lake))+(1:h)/frequency(Lake))
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 1972 9.369630 8.499292 10.239968 8.038563 10.70070
## 1973 8.857415 7.595493 10.119337 6.927472 10.78736
## 1974 8.507123 7.105238 9.909007 6.363126 10.65112
## 1975 8.265032 6.804977 9.725086 6.032071 10.49799
## 1976 8.095244 6.609951 9.580538 5.823684 10.36680
## 1977 7.973771 6.477345 9.470197 5.685185 10.26236
## 1978 7.884584 6.383214 9.385954 5.588436 10.18073
## 1979 7.816971 6.313398 9.320543 5.517454 10.11649
## 1980 7.763773 6.259218 9.268329 5.462755 10.06479
## 1981 7.720210 6.215216 9.225203 5.418520 10.02190
source('https://nmimoto.github.io/R/TS-00.txt')
= Rolling1step.forecast(Lake,
Pred window.size=60,
Arima.order=c(1,0,1),
include.mean=TRUE, # need this for intercept
xreg = TRUE)
##
## Total length 97 , window size 60 .
## Last 37 obs retrospectively forecasted with Rolling 1-step
## prediction using same order and fized window size.
##
## Average Prediction Error: 0.2225
## root Mean Squared Error of Prediction: 0.7817
# Pred
= Arima(Lake, order=c(1,1,2), include.drift=TRUE)
Fit3 Fit3
## Series: Lake
## ARIMA(1,1,2) with drift
##
## Coefficients:
## ar1 ma1 ma2 drift
## 0.6911 -0.6271 -0.3729 -0.0241
## s.e. 0.0956 0.1251 0.1157 0.0100
##
## sigma^2 = 0.4663: log likelihood = -99
## AIC=208.01 AICc=208.67 BIC=220.83
plot(forecast(Fit3, 10))
forecast(Fit3, 10)
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 1972 9.406472 8.527087 10.285857 8.061569 10.75137
## 1973 8.919379 7.629143 10.209615 6.946133 10.89262
## 1974 8.575314 7.124264 10.026365 6.356124 10.79450
## 1975 8.330089 6.804788 9.855390 5.997342 10.66284
## 1976 8.153168 6.591336 9.714999 5.764552 10.54178
## 1977 8.023447 6.442886 9.604009 5.606187 10.44071
## 1978 7.926345 6.335830 9.516860 5.493862 10.35883
## 1979 7.851783 6.255788 9.447778 5.410920 10.29265
## 1980 7.792799 6.193677 9.391920 5.347154 10.23844
## 1981 7.744578 6.143610 9.345547 5.296109 10.19305
source('https://nmimoto.github.io/R/TS-00.txt')
= Rolling1step.forecast(Lake,
Pred window.size=60,
Arima.order=c(1,1,2),
include.drift=TRUE)
##
## Total length 97 , window size 60 .
## Last 37 obs retrospectively forecasted with Rolling 1-step
## prediction using same order and fized window size.
##
## Average Prediction Error: 0.2124
## root Mean Squared Error of Prediction: 0.7781
# Pred
= Arima(Lake, order=c(0,1,0))
Fit4 Fit4
## Series: Lake
## ARIMA(0,1,0)
##
## sigma^2 = 0.5383: log likelihood = -106.49
## AIC=214.98 AICc=215.02 BIC=217.54
plot(forecast(Fit4, 10))
forecast(Fit4, 10)
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 1972 9.96 9.019758 10.90024 8.522024 11.39798
## 1973 9.96 8.630297 11.28970 7.926395 11.99361
## 1974 9.96 8.331453 11.58855 7.469353 12.45065
## 1975 9.96 8.079516 11.84048 7.084048 12.83595
## 1976 9.96 7.857555 12.06244 6.744588 13.17541
## 1977 9.96 7.656887 12.26311 6.437693 13.48231
## 1978 9.96 7.472354 12.44765 6.155473 13.76453
## 1979 9.96 7.300594 12.61941 5.892790 14.02721
## 1980 9.96 7.139274 12.78073 5.646072 14.27393
## 1981 9.96 6.986694 12.93331 5.412721 14.50728
source('https://nmimoto.github.io/R/TS-00.txt')
= Rolling1step.forecast(Lake,
Pred window.size=60,
Arima.order=c(0,1,0),
include.drift=FALSE)
##
## Total length 97 , window size 60 .
## Last 37 obs retrospectively forecasted with Rolling 1-step
## prediction using same order and fized window size.
##
## Average Prediction Error: 0.0843
## root Mean Squared Error of Prediction: 0.7908
# Pred
# From Model1, 2, 3, estimate for $\sigma^2$
= c(0.4784, 0.4612, 0.4663, 0.5383)
Sigma_sq
# Estimates of sigma
sqrt(Sigma_sq)
## [1] 0.6916647 0.6791171 0.6828616 0.7336893
# 0.6916 0.6791 0.6828 0.7336
# Rolling 1-step prediction rMSE for the three models were
# 0.7303, 0.7817, 0.7781, 0.7908