Model pre spanielsko.txt (slajd 14)

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)

Unit root testy

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).

SARIMA model

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

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