library(urca) # unit root testy
library(astsa) # odhadovanie sarima modelov
x <- read.table("http://www.iam.fmph.uniba.sk/institute/stehlikova/cr15/data/spanielsko.txt")
x <- ts(x, frequency=12, start=c(1970,1))
plot(x)
plot(diff(x,12)) # sezonne diferencie
Dve verzie unit root testu: drift
a none
:
summary(ur.df(diff(x,12), type="drift", lags=12, selectlags="BIC"))
##
## ###############################################
## # Augmented Dickey-Fuller Test Unit Root Test #
## ###############################################
##
## Test regression drift
##
##
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1914642 -182358 -21191 137648 1703415
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.408e+04 2.922e+04 2.535 0.012 *
## z.lag.1 -5.490e-01 7.763e-02 -7.073 2.41e-11 ***
## z.diff.lag -1.137e-01 6.980e-02 -1.629 0.105
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 391700 on 203 degrees of freedom
## Multiple R-squared: 0.3185, Adjusted R-squared: 0.3117
## F-statistic: 47.43 on 2 and 203 DF, p-value: < 2.2e-16
##
##
## Value of test-statistic is: -7.0727 25.0115
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau2 -3.46 -2.88 -2.57
## phi1 6.52 4.63 3.81
summary(ur.df(diff(x,12), type="none", lags=12, selectlags="BIC"))
##
## ###############################################
## # Augmented Dickey-Fuller Test Unit Root Test #
## ###############################################
##
## Test regression none
##
##
## Call:
## lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1863075 -94727 40717 209044 1776617
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## z.lag.1 -0.47874 0.07347 -6.517 5.52e-10 ***
## z.diff.lag -0.14916 0.06928 -2.153 0.0325 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 396900 on 204 degrees of freedom
## Multiple R-squared: 0.2969, Adjusted R-squared: 0.29
## F-statistic: 43.07 on 2 and 204 DF, p-value: 2.493e-16
##
##
## Value of test-statistic is: -6.5166
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau1 -2.58 -1.95 -1.62
Modelujeme teda diff(x,12)
.
ACF a PACF:
acf2(diff(x,12))
## ACF PACF
## [1,] 0.37 0.37
## [2,] 0.24 0.12
## [3,] 0.17 0.06
## [4,] 0.15 0.05
## [5,] 0.13 0.04
## [6,] 0.11 0.03
## [7,] 0.12 0.05
## [8,] 0.18 0.11
## [9,] 0.15 0.04
## [10,] 0.10 -0.01
## [11,] 0.01 -0.08
## [12,] -0.25 -0.34
## [13,] 0.03 0.23
## [14,] -0.02 -0.04
## [15,] 0.05 0.09
## [16,] 0.00 -0.04
## [17,] -0.01 -0.04
## [18,] -0.04 -0.05
## [19,] -0.03 0.02
## [20,] -0.10 -0.01
## [21,] -0.12 -0.06
## [22,] -0.10 0.00
## [23,] -0.09 -0.09
## [24,] -0.17 -0.28
## [25,] -0.15 0.15
Pokusy:
# modely tvaru sarima(x,p,0,q,P,1,Q,12)
# co vyslo dobre:
sarima(x,1,0,1,0,1,1,12, details=FALSE)
## $fit
##
## Call:
## stats::arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D,
## Q), period = S), xreg = constant, optim.control = list(trace = trc, REPORT = 1,
## reltol = tol))
##
## Coefficients:
## ar1 ma1 sma1 constant
## 0.7387 -0.3320 -0.6922 10126.511
## s.e. 0.1012 0.1505 0.0602 1673.303
##
## sigma^2 estimated as 1.093e+11: log likelihood = -2922.42, aic = 5854.84
##
## $AIC
## [1] 26.4519
##
## $AICc
## [1] 26.46171
##
## $BIC
## [1] 25.51151
Predikcie na 2 roky (teda 24 mesiacov) z modelu SARIMA(1,0,1)x(0,1,1)_12
sarima.for(x,24,1,0,1,0,1,1,12)
## $pred
## Jan Feb Mar Apr May Jun Jul Aug
## 1989 3613592 4030509 4721301 8178496 9363694
## 1990 2658759 2606071 3093165 3660354 4096803 4802024 8269877 9462949
## 1991 2775380 2723971 3212010
## Sep Oct Nov Dec
## 1989 5633015 3963957 2732035 3449984
## 1990 5738087 4073326 2844578 3564872
## 1991
##
## $se
## Jan Feb Mar Apr May Jun Jul
## 1989 330590.1 356891.1 370457.3 377655.1
## 1990 385902.9 386003.5 386058.3 400499.5 403013.7 404379.1 405122.3
## 1991 405989.8 406000.4 406006.3
## Aug Sep Oct Nov Dec
## 1989 381525.8 383621.6 384760.6 385380.7 385718.6
## 1990 405527.3 405748.2 405868.6 405934.4 405970.2
## 1991