getwd()
## [1] "C:/Users/dell/Desktop/Teaching"
rm(list=ls())
gc()
## used (Mb) gc trigger (Mb) max used (Mb)
## Ncells 293778 7.9 592000 15.9 350000 9.4
## Vcells 323934 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(AirPassengers, package="datasets")
AirPassengers <- as.data.frame(AirPassengers)
ts.plot(AirPassengers$x)
AirPassengers$x
## Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 1949 112 118 132 129 121 135 148 148 136 119 104 118
## 1950 115 126 141 135 125 149 170 170 158 133 114 140
## 1951 145 150 178 163 172 178 199 199 184 162 146 166
## 1952 171 180 193 181 183 218 230 242 209 191 172 194
## 1953 196 196 236 235 229 243 264 272 237 211 180 201
## 1954 204 188 235 227 234 264 302 293 259 229 203 229
## 1955 242 233 267 269 270 315 364 347 312 274 237 278
## 1956 284 277 317 313 318 374 413 405 355 306 271 306
## 1957 315 301 356 348 355 422 465 467 404 347 305 336
## 1958 340 318 362 348 363 435 491 505 404 359 310 337
## 1959 360 342 406 396 420 472 548 559 463 407 362 405
## 1960 417 391 419 461 472 535 622 606 508 461 390 432
bulkfit(AirPassengers$x)
## $res
## ar d ma AIC
## [1,] 0 0 0 1790.368
## [2,] 0 0 1 1618.863
## [3,] 0 0 2 1522.122
## [4,] 0 1 0 1413.909
## [5,] 0 1 1 1397.258
## [6,] 0 1 2 1397.093
## [7,] 0 2 0 1450.596
## [8,] 0 2 1 1411.368
## [9,] 0 2 2 1394.373
## [10,] 1 0 0 1428.179
## [11,] 1 0 1 1409.748
## [12,] 1 0 2 1411.050
## [13,] 1 1 0 1401.853
## [14,] 1 1 1 1394.683
## [15,] 1 1 2 1385.497
## [16,] 1 2 0 1447.028
## [17,] 1 2 1 1398.929
## [18,] 1 2 2 1391.910
## [19,] 2 0 0 1413.639
## [20,] 2 0 1 1408.249
## [21,] 2 0 2 1408.343
## [22,] 2 1 0 1396.588
## [23,] 2 1 1 1378.338
## [24,] 2 1 2 1387.409
## [25,] 2 2 0 1440.078
## [26,] 2 2 1 1393.882
## [27,] 2 2 2 1392.659
##
## $min
## ar d ma AIC
## 2.000 1.000 1.000 1378.338
ArimaModel.1 <-
Arima(AirPassengers$x,order=c(2,
1,1),
include.mean=1,
seasonal=list(order=c(0,0,0),
period=12))
ArimaModel.1
## Series: AirPassengers$x
## ARIMA(2,1,1)
##
## Coefficients:
## ar1 ar2 ma1
## 1.0906 -0.4890 -0.8438
## s.e. 0.0776 0.0744 0.0427
##
## sigma^2 estimated as 844.6: log likelihood=-685.17
## AIC=1378.34 AICc=1378.63 BIC=1390.19
ArimaModel.2 <-
Arima(AirPassengers$x,order=c(2,
1,1),
include.mean=0,
seasonal=list(order=c(0,0,0),
period=12))
ArimaModel.2
## Series: AirPassengers$x
## ARIMA(2,1,1)
##
## Coefficients:
## ar1 ar2 ma1
## 1.0906 -0.4890 -0.8438
## s.e. 0.0776 0.0744 0.0427
##
## sigma^2 estimated as 844.6: log likelihood=-685.17
## AIC=1378.34 AICc=1378.63 BIC=1390.19
HoltWintersModel.3 <-
HoltWinters(AirPassengers$x,
seasonal='add')
HoltWintersModel.3
## Holt-Winters exponential smoothing with trend and additive seasonal component.
##
## Call:
## HoltWinters(x = AirPassengers$x, seasonal = "add")
##
## Smoothing parameters:
## alpha: 0.2479595
## beta : 0.03453373
## gamma: 1
##
## Coefficients:
## [,1]
## a 477.827781
## b 3.127627
## s1 -27.457685
## s2 -54.692464
## s3 -20.174608
## s4 12.919120
## s5 18.873607
## s6 75.294426
## s7 152.888368
## s8 134.613464
## s9 33.778349
## s10 -18.379060
## s11 -87.772408
## s12 -45.827781
HoltWintersModel.4 <-
HoltWinters(AirPassengers$x,
gamma=0,beta=0)
HoltWintersModel.4
## Holt-Winters exponential smoothing with trend and additive seasonal component.
##
## Call:
## HoltWinters(x = AirPassengers$x, beta = 0, gamma = 0)
##
## Smoothing parameters:
## alpha: 0.9999339
## beta : 0
## gamma: 0
##
## Coefficients:
## [,1]
## a 442.31761986
## b 1.14568765
## s1 -14.81944444
## s2 -5.65277778
## s3 7.51388889
## s4 0.01388889
## s5 -10.98611111
## s6 11.68055556
## s7 22.63888889
## s8 22.18055556
## s9 9.47222222
## s10 -8.15277778
## s11 -23.56944444
## s12 -10.31944444
a=ets(AirPassengers$x)
a
## ETS(M,A,M)
##
## Call:
## ets(y = AirPassengers$x)
##
## Smoothing parameters:
## alpha = 0.5901
## beta = 0.0058
## gamma = 1e-04
##
## Initial states:
## l = 126.9791
## b = 1.6483
## s=0.8865 0.7928 0.9226 1.0582 1.2186 1.2371
## 1.1069 0.9818 0.985 1.0149 0.8946 0.901
##
## sigma: 0.0367
##
## AIC AICc BIC
## 1391.174 1395.457 1438.691
predict(a,10)
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Jan 1961 444.9979 424.0692 465.9267 412.9901 477.0057
## Feb 1961 444.1187 419.8324 468.4049 406.9760 481.2613
## Mar 1961 506.4125 475.2920 537.5330 458.8178 554.0072
## Apr 1961 493.9890 460.5937 527.3844 442.9153 545.0628
## May 1961 494.8520 458.5835 531.1205 439.3841 550.3198
## Jun 1961 560.6983 516.6174 604.7792 493.2823 628.1143
## Jul 1961 629.8161 577.1310 682.5012 549.2412 710.3910
## Aug 1961 623.4715 568.3301 678.6129 539.1400 707.8031
## Sep 1961 544.0719 493.4577 594.6861 466.6642 621.4796
## Oct 1961 476.6947 430.2464 523.1430 405.6581 547.7313
str(predict(a,20))
## List of 9
## $ model :List of 18
## ..$ loglik : num -680
## ..$ aic : num 1391
## ..$ bic : num 1439
## ..$ aicc : num 1395
## ..$ mse : num 121
## ..$ amse : num 194
## ..$ fit :List of 4
## .. ..$ value : num 1359
## .. ..$ par : num [1:16] 5.90e-01 5.82e-03 1.21e-04 1.27e+02 1.65 ...
## .. ..$ fail : int 1
## .. ..$ fncount: int 2001
## ..$ residuals : Time-Series [1:144] from 1949 to 1961: -0.0336 0.0329 -0.0134 -0.011 -0.0747 ...
## ..$ fitted : Time-Series [1:144] from 1949 to 1961: 116 114 134 130 131 ...
## ..$ states : mts [1:145, 1:14] 127 126 130 131 132 ...
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : NULL
## .. .. ..$ : chr [1:14] "l" "b" "s1" "s2" ...
## .. ..- attr(*, "tsp")= num [1:3] 1949 1961 12
## .. ..- attr(*, "class")= chr [1:3] "mts" "ts" "matrix"
## ..$ par : Named num [1:16] 5.90e-01 5.82e-03 1.21e-04 1.27e+02 1.65 ...
## .. ..- attr(*, "names")= chr [1:16] "alpha" "beta" "gamma" "l" ...
## ..$ m : num 12
## ..$ method : chr "ETS(M,A,M)"
## ..$ components: chr [1:4] "M" "A" "M" "FALSE"
## ..$ call : language ets(y = AirPassengers$x)
## ..$ initstate : Named num [1:14] 126.979 1.648 0.887 0.793 0.923 ...
## .. ..- attr(*, "names")= chr [1:14] "l" "b" "s1" "s2" ...
## ..$ sigma2 : num 0.00135
## ..$ x : Time-Series [1:144] from 1949 to 1961: 112 118 132 129 121 135 148 148 136 119 ...
## ..- attr(*, "class")= chr "ets"
## $ mean : Time-Series [1:20] from 1961 to 1963: 445 444 506 494 495 ...
## $ level : num [1:2] 80 95
## $ x : Time-Series [1:144] from 1949 to 1961: 112 118 132 129 121 135 148 148 136 119 ...
## $ upper : mts [1:20, 1:2] 466 468 538 527 531 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : NULL
## .. ..$ : chr [1:2] "Series 1" "Series 2"
## ..- attr(*, "tsp")= num [1:3] 1961 1963 12
## ..- attr(*, "class")= chr [1:3] "mts" "ts" "matrix"
## $ lower : mts [1:20, 1:2] 424 420 475 461 459 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : NULL
## .. ..$ : chr [1:2] "Series 1" "Series 2"
## ..- attr(*, "tsp")= num [1:3] 1961 1963 12
## ..- attr(*, "class")= chr [1:3] "mts" "ts" "matrix"
## $ fitted : Time-Series [1:144] from 1949 to 1961: 116 114 134 130 131 ...
## $ method : chr "ETS(M,A,M)"
## $ residuals: Time-Series [1:144] from 1949 to 1961: -0.0336 0.0329 -0.0134 -0.011 -0.0747 ...
## - attr(*, "class")= chr "forecast"
b=auto.arima(AirPassengers$x)
b
## Series: AirPassengers$x
## ARIMA(0,1,1)(0,1,0)[12]
##
## Coefficients:
## ma1
## -0.3184
## s.e. 0.0877
##
## sigma^2 estimated as 137.3: log likelihood=-508.32
## AIC=1020.64 AICc=1020.73 BIC=1026.39
predict(b,10)
## $pred
## Jan Feb Mar Apr May Jun Jul
## 1961 446.7582 420.7582 448.7582 490.7582 501.7582 564.7582 651.7582
## Aug Sep Oct
## 1961 635.7582 537.7582 490.7582
##
## $se
## Jan Feb Mar Apr May Jun Jul
## 1961 11.71600 14.17842 16.27239 18.12605 19.80699 21.35602 22.80006
## Aug Sep Oct
## 1961 24.15793 25.44344 26.66705
newdataset=predar3(b,fore1=48)
newdataset2=predar3(ArimaModel.1,fore1=48)
