getwd()
## [1] "C:/Users/dell/Desktop/Teaching"
rm(list=ls())
gc()
## used (Mb) gc trigger (Mb) max used (Mb)
## Ncells 293614 7.9 460000 12.3 350000 9.4
## Vcells 323050 2.5 786432 6.0 677529 5.2
sessionInfo()
## R version 3.2.2 (2015-08-14)
## Platform: i386-w64-mingw32/i386 (32-bit)
## Running under: Windows 7 (build 7601) Service Pack 1
##
## locale:
## [1] LC_COLLATE=English_United States.1252
## [2] LC_CTYPE=English_United States.1252
## [3] LC_MONETARY=English_United States.1252
## [4] LC_NUMERIC=C
## [5] LC_TIME=English_United States.1252
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## loaded via a namespace (and not attached):
## [1] magrittr_1.5 tools_3.2.2 htmltools_0.3 stringi_1.0-1 rmarkdown_0.7
## [6] knitr_1.10.5 stringr_1.0.0 digest_0.6.8 evaluate_0.7
memory.limit()
## [1] 1535
library(Rcmdr)
## Loading required package: splines
## Loading required package: RcmdrMisc
## Loading required package: car
## Loading required package: sandwich
## The Commander GUI is launched only in interactive sessions
library(RcmdrPlugin.epack)
## Loading required package: TeachingDemos
## Loading required package: tseries
## Loading required package: abind
## Loading required package: MASS
## 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: forecast
## Loading required package: timeDate
## This is forecast 6.2
data(woolyrnq, package="forecast")
woolyrnq <- as.data.frame(woolyrnq)
a=ets(woolyrnq$x)
a
## ETS(M,N,A)
##
## Call:
## ets(y = woolyrnq$x)
##
## Smoothing parameters:
## alpha = 0.9374
## gamma = 1e-04
##
## Initial states:
## l = 6765.8183
## s=27.5251 463.9382 128.9699 -620.4332
##
## sigma: 0.0748
##
## AIC AICc BIC
## 2016.427 2017.177 2033.102
predict(a,10)
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 1994 Q4 5962.242 5390.649 6533.835 5088.066 6836.418
## 1995 Q1 5314.252 4573.814 6054.690 4181.849 6446.654
## 1995 Q2 6063.343 5137.379 6989.307 4647.203 7479.482
## 1995 Q3 6398.296 5304.278 7492.314 4725.139 8071.452
## 1995 Q4 5962.242 4744.109 7180.375 4099.269 7825.215
## 1996 Q1 5314.252 4006.171 6622.332 3313.715 7314.788
## 1996 Q2 6063.343 4639.986 7486.699 3886.507 8240.179
## 1996 Q3 6398.296 4858.444 7938.148 4043.296 8753.296
## 1996 Q4 5962.242 4330.095 7594.389 3466.088 8458.396
## 1997 Q1 5314.252 3612.211 7016.292 2711.205 7917.298
str(predict(a,20))
## List of 9
## $ model :List of 18
## ..$ loglik : num -1002
## ..$ aic : num 2016
## ..$ bic : num 2033
## ..$ aicc : num 2017
## ..$ mse : num 195522
## ..$ amse : num 310159
## ..$ fit :List of 4
## .. ..$ value : num 2004
## .. ..$ par : num [1:6] 9.37e-01 1.00e-04 6.77e+03 2.75e+01 4.64e+02 ...
## .. ..$ fail : int 0
## .. ..$ fncount: int 689
## ..$ residuals : Time-Series [1:119] from 1965 to 1994: 0.00433 -0.03045 -0.0601 0.0702 0.13389 ...
## ..$ fitted : Time-Series [1:119] from 1965 to 1994: 6145 6920 7057 6223 5985 ...
## ..$ states : mts [1:120, 1:5] 6766 6791 6593 6196 6605 ...
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : NULL
## .. .. ..$ : chr [1:5] "l" "s1" "s2" "s3" ...
## .. ..- attr(*, "tsp")= num [1:3] 1965 1994 4
## .. ..- attr(*, "class")= chr [1:3] "mts" "ts" "matrix"
## ..$ par : Named num [1:6] 9.37e-01 1.00e-04 6.77e+03 2.75e+01 4.64e+02 ...
## .. ..- attr(*, "names")= chr [1:6] "alpha" "gamma" "l" "s0" ...
## ..$ m : num 4
## ..$ method : chr "ETS(M,N,A)"
## ..$ components: chr [1:4] "M" "N" "A" "FALSE"
## ..$ call : language ets(y = woolyrnq$x)
## ..$ initstate : Named num [1:5] 6765.8 27.5 463.9 129 -620.4
## .. ..- attr(*, "names")= chr [1:5] "l" "s1" "s2" "s3" ...
## ..$ sigma2 : num 0.0056
## ..$ x : Time-Series [1:119] from 1965 to 1994: 6172 6709 6633 6660 6786 ...
## ..- attr(*, "class")= chr "ets"
## $ mean : Time-Series [1:20] from 1995 to 2000: 5962 5314 6063 6398 5962 ...
## $ level : num [1:2] 80 95
## $ x : Time-Series [1:119] from 1965 to 1994: 6172 6709 6633 6660 6786 ...
## $ upper : mts [1:20, 1:2] 6534 6055 6989 7492 7180 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : NULL
## .. ..$ : chr [1:2] "Series 1" "Series 2"
## ..- attr(*, "tsp")= num [1:3] 1995 2000 4
## ..- attr(*, "class")= chr [1:3] "mts" "ts" "matrix"
## $ lower : mts [1:20, 1:2] 5391 4574 5137 5304 4744 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : NULL
## .. ..$ : chr [1:2] "Series 1" "Series 2"
## ..- attr(*, "tsp")= num [1:3] 1995 2000 4
## ..- attr(*, "class")= chr [1:3] "mts" "ts" "matrix"
## $ fitted : Time-Series [1:119] from 1965 to 1994: 6145 6920 7057 6223 5985 ...
## $ method : chr "ETS(M,N,A)"
## $ residuals: Time-Series [1:119] from 1965 to 1994: 0.00433 -0.03045 -0.0601 0.0702 0.13389 ...
## - attr(*, "class")= chr "forecast"
b=auto.arima(woolyrnq$x)
b
## Series: woolyrnq$x
## ARIMA(1,0,0)(0,1,1)[4]
##
## Coefficients:
## ar1 sma1
## 0.8077 -0.6669
## s.e. 0.0629 0.0944
##
## sigma^2 estimated as 172819: log likelihood=-858
## AIC=1722 AICc=1722.21 BIC=1730.23
predict(b,10)
## $pred
## Qtr1 Qtr2 Qtr3 Qtr4
## 1994 5798.859
## 1995 5046.575 5843.460 6132.798 5586.269
## 1996 4874.866 5704.769 6020.777 5495.789
## 1997 4801.785
##
## $se
## Qtr1 Qtr2 Qtr3 Qtr4
## 1994 415.7154
## 1995 534.3820 599.2640 638.0458 711.7536
## 1996 755.9759 783.4817 800.9171 846.0810
## 1997 874.2886
newdataset=predar3(b,fore1=48)
bulkfit(woolyrnq$x)
## $res
## ar d ma AIC
## [1,] 0 0 0 1994.256
## [2,] 0 0 1 1918.450
## [3,] 0 0 2 1878.005
## [4,] 0 1 0 1875.308
## [5,] 0 1 1 1860.141
## [6,] 0 1 2 1842.926
## [7,] 0 2 0 1950.905
## [8,] 0 2 1 1867.198
## [9,] 0 2 2 1853.993
## [10,] 1 0 0 1882.118
## [11,] 1 0 1 1882.007
## [12,] 1 0 2 1879.646
## [13,] 1 1 0 1876.419
## [14,] 1 1 1 1855.804
## [15,] 1 1 2 1844.874
## [16,] 1 2 0 1944.290
## [17,] 1 2 1 1868.474
## [18,] 1 2 2 1869.716
## [19,] 2 0 0 1884.058
## [20,] 2 0 1 1882.261
## [21,] 2 0 2 1875.342
## [22,] 2 1 0 1821.890
## [23,] 2 1 1 1816.072
## [24,] 2 1 2 1784.722
## [25,] 2 2 0 1897.500
## [26,] 2 2 1 1815.369
## [27,] 2 2 2 1810.019
##
## $min
## ar d ma AIC
## 2.000 1.000 2.000 1784.722
ArimaModel.5 <- Arima(woolyrnq$x,
order=c(2,1,2),
include.mean=1,
seasonal=list(order=c(0,0,0),
period=4))
ArimaModel.5
## Series: woolyrnq$x
## ARIMA(2,1,2)
##
## Coefficients:
## ar1 ar2 ma1 ma2
## 0.0076 -0.9791 -0.0699 0.7808
## s.e. 0.0197 0.0184 0.0730 0.0738
##
## sigma^2 estimated as 194860: log likelihood=-887.36
## AIC=1784.72 AICc=1785.26 BIC=1798.58
newdataset2=predar3(ArimaModel.5,fore1=48)
