Class Web Page



Lake Huron - Direct fit with ARMA

  # Direct fit of ARMA(p,q) to Lake Huron data
  library(forecast)
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
  source('https://nmimoto.github.io/R/TS-00.txt')
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
  Huron = read.csv("https://nmimoto.github.io/datasets/lake.csv")
  Y = ts(Huron[,1], start=c(1875,1), freq=1)
  plot(Y)

  layout(matrix(1:2, 1,2))
  acf(Y)
  pacf(Y)

  layout(1)

  Fit1 = auto.arima(Y, d=0, stepwise=FALSE, approximation=FALSE)
  Fit1              # best ARMA(p,q) by AICc
## Series: Y 
## 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)
## 
##     'tseries' version: 0.10-49
## 
##     'tseries' is a package for time series analysis and
##     computational finance.
## 
##     See 'library(help="tseries")' for details.

##   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

ARMA(1,1) with mean is the final model. \[\begin{align*} Y_t \mbox{ is the observation.} \\\\ Y_t = \mu + X_t \\\\ X_t = \phi_1 X_{t-1} + \epsilon_t + \theta_1 \epsilon_{t-1} \end{align*}\] Parameter estimates are: \[\begin{align*} \hat \mu &= 9.1290 \\ \hat \phi_1 &= 0.7665 \\ \hat \theta_1 &= 0.3393 \\ \hat \sigma^2 &= \sqrt{0.4784} = 0.6917 \end{align*}\]

Static 10-step Forecast

  # Using Fit1 from above
  plot(forecast(Fit1, 10))
  abline(h=Fit1$coef["intercept"])

  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



(Retroactive) Rolling 1-step Forecast

  # Using Y from above
  source('https://nmimoto.github.io/R/TS-00.txt')

  Pred = Rolling1step.forecast(Y,
                               window.size=70,
                               Arima.order=c(1,0,1),
                               include.mean=TRUE)
## 
##   Total length 97 , window size 70 .
##   Last 27 obs retrospectively forecasted with Rolling 1-step
##       prediction using same order and fized window size.
## 
##   Average Prediction Error:  -0.0512
##   root Mean Squared Error of Prediction:   0.771

  Pred
## $Actual
## Time Series:
## Start = 1945 
## End = 1971 
## Frequency = 1 
##  [1]  9.22  9.38  9.10  7.95  8.12  9.75 10.85 10.41  9.96  9.61
## [11]  8.76  8.18  7.21  7.13  9.10  8.25  7.91  6.89  5.96  6.80
## [21]  7.68  8.38  8.52  9.74  9.31  9.89  9.96
## 
## $Predicted
## Time Series:
## Start = 1945 
## End = 1971 
## Frequency = 1 
##  [1]  9.413844  9.125644  9.426446  8.956047  7.789492  8.490787
##  [7] 10.123483 10.707349  9.881199  9.714921  9.364789  8.523881
## [13]  8.272220  7.285383  7.705719  9.700834  7.711935  8.375328
## [19]  7.000006  6.533042  7.626733  8.064683  8.618905  8.487383
## [25]  9.906468  8.717502  9.877759