Check the help files of functions
?acf
?pacf
?astsa::sarima
?astsa::sarima.for
?astsa::cmort
f = astsa::sarima(astsa::cmort, p=1, d=0, q=0)
## initial value 2.301627
## iter 2 value 1.847932
## iter 3 value 1.847931
## iter 4 value 1.847930
## iter 5 value 1.847930
## iter 6 value 1.847928
## iter 7 value 1.847927
## iter 8 value 1.847927
## iter 9 value 1.847927
## iter 10 value 1.847926
## iter 11 value 1.847925
## iter 12 value 1.847925
## iter 13 value 1.847925
## iter 13 value 1.847925
## iter 13 value 1.847925
## final value 1.847925
## converged
## initial value 1.848674
## iter 2 value 1.848674
## iter 3 value 1.848668
## iter 4 value 1.848666
## iter 5 value 1.848665
## iter 6 value 1.848664
## iter 6 value 1.848664
## final value 1.848664
## converged
str(f)
## List of 6
## $ fit :List of 14
## ..$ coef : Named num [1:2] 0.771 88.739
## .. ..- attr(*, "names")= chr [1:2] "ar1" "xmean"
## ..$ sigma2 : num 40.3
## ..$ var.coef : num [1:2, 1:2] 0.000789 0.000181 0.000181 1.497858
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : chr [1:2] "ar1" "xmean"
## .. .. ..$ : chr [1:2] "ar1" "xmean"
## ..$ mask : logi [1:2] TRUE TRUE
## ..$ loglik : num -1660
## ..$ aic : num 3326
## ..$ arma : int [1:7] 1 0 0 0 -1 0 0
## ..$ residuals: Time-Series [1:508] from 1970 to 1980: 5.797 8.8721 -6.6461 4.9746 -0.0722 ...
## ..$ call : language stats::arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D, Q), period = S), xreg = xmean, i| __truncated__ ...
## ..$ series : chr "xdata"
## ..$ code : int 0
## ..$ n.cond : int 0
## ..$ nobs : int 508
## ..$ model :List of 10
## .. ..$ phi : num 0.771
## .. ..$ theta: num(0)
## .. ..$ Delta: num(0)
## .. ..$ Z : num 1
## .. ..$ a : num -3.25
## .. ..$ P : num [1, 1] 0
## .. ..$ T : num [1, 1] 0.771
## .. ..$ V : num [1, 1] 1
## .. ..$ h : num 0
## .. ..$ Pn : num [1, 1] 1
## ..- attr(*, "class")= chr "Arima"
## $ degrees_of_freedom: int 506
## $ ttable : num [1:2, 1:4] 0.7715 88.7392 0.0281 1.2239 27.4731 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:2] "ar1" "xmean"
## .. ..$ : chr [1:4] "Estimate" "SE" "t.value" "p.value"
## $ AIC : num 6.55
## $ AICc : num 6.55
## $ BIC : num 6.57
f$ttable
## Estimate SE t.value p.value
## ar1 0.7715 0.0281 27.4731 0
## xmean 88.7392 1.2239 72.5070 0
p = astsa::sarima.for(astsa::cmort, n.ahead=3, p=1, d=0, q=0)
str(p)
## List of 2
## $ pred: Time-Series [1:3] from 1980 to 1980: 86.2 86.8 87.2
## $ se : Time-Series [1:3] from 1980 to 1980: 6.35 8.01 8.86
df = read.csv("HW5_Q2_data.csv")
x = ts(df$x)
str(x)
## Time-Series [1:200] from 1 to 200: 0.358 -1.965 -3.7 -3.997 -2.596 ...
astsa::sarima(x, p=1, d=0, q=1)
## initial value 0.644518
## iter 2 value 0.297716
## iter 3 value 0.069358
## iter 4 value 0.064231
## iter 5 value 0.053791
## iter 6 value 0.048103
## iter 7 value 0.048049
## iter 8 value 0.048038
## iter 9 value 0.048033
## iter 10 value 0.048014
## iter 11 value 0.048001
## iter 12 value 0.048001
## iter 13 value 0.048001
## iter 13 value 0.048001
## iter 13 value 0.048001
## final value 0.048001
## converged
## initial value 0.048408
## iter 2 value 0.048398
## iter 3 value 0.048388
## iter 4 value 0.048380
## iter 5 value 0.048377
## iter 6 value 0.048373
## iter 7 value 0.048372
## iter 8 value 0.048371
## iter 9 value 0.048371
## iter 10 value 0.048371
## iter 11 value 0.048371
## iter 12 value 0.048371
## iter 12 value 0.048371
## iter 12 value 0.048371
## final value 0.048371
## converged
## $fit
##
## Call:
## stats::arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D,
## Q), period = S), xreg = xmean, include.mean = FALSE, transform.pars = trans,
## fixed = fixed, optim.control = list(trace = trc, REPORT = 1, reltol = tol))
##
## Coefficients:
## ar1 ma1 xmean
## 0.7062 0.3850 0.1235
## s.e. 0.0606 0.0809 0.3443
##
## sigma^2 estimated as 1.094: log likelihood = -293.46, aic = 594.92
##
## $degrees_of_freedom
## [1] 197
##
## $ttable
## Estimate SE t.value p.value
## ar1 0.7062 0.0606 11.6553 0.0000
## ma1 0.3850 0.0809 4.7592 0.0000
## xmean 0.1235 0.3443 0.3586 0.7203
##
## $AIC
## [1] 2.97462
##
## $AICc
## [1] 2.975232
##
## $BIC
## [1] 3.040586
arima(x, order = c(1, 0, 1), fixed = c(0.7, 0.4, 0.0))
##
## Call:
## arima(x = x, order = c(1, 0, 1), fixed = c(0.7, 0.4, 0))
##
## Coefficients:
## ar1 ma1 intercept
## 0.7 0.4 0.0
##
## sigma^2 estimated as 1.095: log likelihood = -293.55, aic = 589.09
astsa::sarima.for(x, n.ahead=3, p=1, d=0, q=1, fixed = c(0.7, 0.4, 0.0))
## $pred
## Time Series:
## Start = 201
## End = 203
## Frequency = 1
## [1] -3.429269 -2.400488 -1.680342
##
## $se
## Time Series:
## Start = 201
## End = 203
## Frequency = 1
## [1] 1.046489 1.555718 1.752018
If the true model is \[ \begin{aligned} X_{t+1} &= \varphi X_t + W_{t+1} \\ X_{t+2} &= \varphi X_{t+1} + W_{t+2} \end{aligned} \] then \[ \begin{aligned} \text{E}[X_{t+1}| X_t ] &= \text{E}[\varphi X_t + W_{t+1} | X_t] \\ &= \varphi \text{E}[X_t| X_t] + \text{E}[W_{t+1} | X_t] \\ &= \varphi X_t \end{aligned} \] and also \[ \begin{aligned} \text{E}[X_{t+2}| X_t ] &= \text{E}[ \varphi X_{t+1} + W_{t+2}| X_t ] \\ &= \varphi \text{E}[ X_{t+1}| X_t ] + \text{E}[W_{t+2}| X_t ] \\ &= \varphi^2 X_t \end{aligned} \]
See! Math is easy!