1
library(astsa)
y <- window(ts(econ5[, 1]), start = 81, end = 150)
plot(y)
# linearny trend -> diferencujeme
# pozrieme si diferencie
plot(diff(y))
abline(h=0, col = "red")
abline(h=mean(diff(y)), col = "blue")
# aby povodne data mali linearny trend, v diferenciach treba konstantny clen
# teda diff(y) maju nenulovu strednu hodnotu
# treba zistit, ci je v nich jednotkovy koren
# -> ADF test
# diff(y) s nenulovou strednou hodnotou
# bez trendu
# -> type = "drift"
# skusme lags=5 (maximalny pocet ktory ma R uvazovat)
# BIC je standard
library(urca)
ur.df(diff(y), type = "drift", lags = 5, selectlags = "BIC")
##
## ###############################################################
## # Augmented Dickey-Fuller Test Unit Root / Cointegration Test #
## ###############################################################
##
## The value of the test statistic is: -3.5473 6.2948
# toto da len statistiku -> treba summary
summary(ur.df(diff(y), type = "drift", lags = 5, 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
## -0.80333 -0.22357 0.01922 0.22106 1.80569
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.02351 0.05155 0.456 0.650038
## z.lag.1 -0.44570 0.12565 -3.547 0.000763 ***
## z.diff.lag -0.06349 0.12886 -0.493 0.624015
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4045 on 60 degrees of freedom
## Multiple R-squared: 0.2413, Adjusted R-squared: 0.2161
## F-statistic: 9.544 on 2 and 60 DF, p-value: 0.000252
##
##
## Value of test-statistic is: -3.5473 6.2948
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau2 -3.51 -2.89 -2.58
## phi1 6.70 4.71 3.86
# statistika -3.547
# kriticka hodnota -2.89
# statistika < kriticka hodnota -> H0 sa zamieta -> jednotkovy koren sa zamieta
# H1 je stacionarita -> data su stacionarne
# -> ARMA model pre diff(y)
# sarima(y, p, 1, q)
library(astsa)
acf2(diff(y)) # naraz ACF aj PACF
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
## ACF 0.52 0.32 0.08 -0.13 -0.11 -0.20 -0.14 -0.21 -0.18 -0.10 -0.14 -0.18
## PACF 0.52 0.07 -0.15 -0.19 0.07 -0.13 -0.01 -0.17 -0.01 0.02 -0.12 -0.22
## [,13] [,14] [,15] [,16] [,17] [,18] [,19]
## ACF -0.18 -0.16 -0.13 -0.03 0.00 0.02 0.10
## PACF -0.04 -0.04 -0.10 -0.02 -0.06 -0.08 0.06
sarima(y, 1, 1, 0) # chceme konstatny clen kvoli trendu v podovnych datach
## initial value -0.799245
## iter 2 value -0.954532
## iter 3 value -0.954534
## iter 4 value -0.954536
## iter 4 value -0.954536
## iter 4 value -0.954536
## final value -0.954536
## converged
## initial value -0.958768
## iter 2 value -0.958822
## iter 3 value -0.958834
## iter 4 value -0.958834
## iter 5 value -0.958834
## iter 5 value -0.958834
## iter 5 value -0.958834
## final value -0.958834
## converged
## $fit
##
## Call:
## arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D, Q), period = S),
## xreg = constant, transform.pars = trans, fixed = fixed, optim.control = list(trace = trc,
## REPORT = 1, reltol = tol))
##
## Coefficients:
## ar1 constant
## 0.5098 0.0564
## s.e. 0.1019 0.0926
##
## sigma^2 estimated as 0.1463: log likelihood = -31.75, aic = 69.49
##
## $degrees_of_freedom
## [1] 67
##
## $ttable
## Estimate SE t.value p.value
## ar1 0.5098 0.1019 5.0038 0.0000
## constant 0.0564 0.0926 0.6092 0.5444
##
## $AIC
## [1] 1.007166
##
## $AICc
## [1] 1.009801
##
## $BIC
## [1] 1.104301
# rezidua OK
# treba overit stacionaritu a invertovatelnost
2
library(astsa)
y <- ts(window(gas, start = c(2006,2), end = c(2009,52)))
plot(y)
summary(ur.df(y, type = "drift", lags = 5, 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
## -25.3462 -4.9850 0.2505 6.2098 20.4904
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.08163 2.55761 2.378 0.0184 *
## z.lag.1 -0.03015 0.01236 -2.439 0.0156 *
## z.diff.lag1 0.19176 0.06971 2.751 0.0065 **
## z.diff.lag2 0.17296 0.07013 2.466 0.0145 *
## z.diff.lag3 0.16363 0.06980 2.344 0.0201 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 9.203 on 196 degrees of freedom
## Multiple R-squared: 0.1452, Adjusted R-squared: 0.1277
## F-statistic: 8.322 on 4 and 196 DF, p-value: 3.21e-06
##
##
## Value of test-statistic is: -2.4392 2.9776
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau2 -3.46 -2.88 -2.57
## phi1 6.52 4.63 3.81
# stat > krit. hodnota -> H0 nezamietame -> nezamietame nulovy koren
# -> data diferencujem
plot(diff(y))
# nulova stredna hodnota
# tvrdili sme ze y nema linearny trend
summary(ur.df(diff(y), type = "none", lags = 5, 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
## -25.220 -5.945 1.077 6.739 21.409
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## z.lag.1 -0.60103 0.08465 -7.100 2.18e-11 ***
## z.diff.lag -0.18336 0.06975 -2.629 0.00924 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 9.373 on 198 degrees of freedom
## Multiple R-squared: 0.3901, Adjusted R-squared: 0.384
## F-statistic: 63.33 on 2 and 198 DF, p-value: < 2.2e-16
##
##
## Value of test-statistic is: -7.1004
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau1 -2.58 -1.95 -1.62
# zamietame jednotkovy koren -> ARMA model pre diferencie
# sarima(y, p, 1, q, no.constant = TRUE)
# model pre diferencie nema mat konstantu -> no.constant=TRUE
acf2(diff(y))
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
## ACF 0.25 0.23 0.23 -0.01 0.07 0.04 -0.01 0.15 -0.02 0.04 0.10 0 -0.02
## PACF 0.25 0.18 0.15 -0.14 0.03 0.02 0.00 0.14 -0.09 0.02 0.06 0 -0.08
## [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24] [,25]
## ACF -0.02 -0.03 0.02 -0.02 0.06 -0.04 -0.03 -0.08 -0.12 -0.24 -0.18 -0.21
## PACF -0.04 0.02 0.03 -0.01 0.06 -0.11 0.00 -0.07 -0.06 -0.22 -0.04 -0.07
sarima(y, 0, 1, 3, no.constant = TRUE)
## initial value 2.283784
## iter 2 value 2.217640
## iter 3 value 2.216156
## iter 4 value 2.215974
## iter 5 value 2.215969
## iter 6 value 2.215969
## iter 6 value 2.215969
## iter 6 value 2.215969
## final value 2.215969
## converged
## initial value 2.216032
## iter 2 value 2.216029
## iter 3 value 2.216029
## iter 3 value 2.216029
## iter 3 value 2.216029
## final value 2.216029
## converged
## $fit
##
## Call:
## arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D, Q), period = S),
## include.mean = !no.constant, transform.pars = trans, fixed = fixed, optim.control = list(trace = trc,
## REPORT = 1, reltol = tol))
##
## Coefficients:
## ma1 ma2 ma3
## 0.2039 0.2128 0.2418
## s.e. 0.0673 0.0696 0.0662
##
## sigma^2 estimated as 84: log likelihood = -748.8, aic = 1505.61
##
## $degrees_of_freedom
## [1] 203
##
## $ttable
## Estimate SE t.value p.value
## ma1 0.2039 0.0673 3.0276 0.0028
## ma2 0.2128 0.0696 3.0555 0.0025
## ma3 0.2418 0.0662 3.6531 0.0003
##
## $AIC
## [1] 7.30877
##
## $AICc
## [1] 7.309347
##
## $BIC
## [1] 7.373389
# rezidua OK
# treba overit stacionaritu a invertovatelnost