#2.Визуальный анализ временного ряда
plot(Sample1$Ряд3, type="l")
# 3.Стационарность=============================================================================================
# 1)Сначала строим график частичной автокорреляции, он нужен для определения порядка тестирования гипотез о стационарности
pacf(Sample1$Ряд3)
# 2)Видим по графику, что значим 1 лаг => в следующей функции будет lags=1
# 3)Тест Дики-Фуллера(type=drift, это значит, что изменения во времени происходят без явного тренда)
df <- ur.df(Sample1$Ряд3, type='drift', lags=1)
summary(df)
##
## ###############################################
## # 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
## -10.5669 -1.8001 0.0114 1.7502 11.7886
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.101875 0.098393 11.199 <2e-16 ***
## z.lag.1 -1.411078 0.054486 -25.898 <2e-16 ***
## z.diff.lag 0.000563 0.032408 0.017 0.986
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.738 on 951 degrees of freedom
## Multiple R-squared: 0.705, Adjusted R-squared: 0.7044
## F-statistic: 1137 on 2 and 951 DF, p-value: < 2.2e-16
##
##
## Value of test-statistic is: -25.898 335.3533
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau2 -3.43 -2.86 -2.57
## phi1 6.43 4.59 3.78
# 4) В выводе смотрим на две вещи: Value of test-statistic и Critical values for test statistics.Если расчетное больше критического (по модулю), то гипотезу Null о нестационарности отвергаем
# В данном случае -25.898=test и -3.43=critical ->
#Null о нестационарности отвергаем -> ряд стационарный, к первым разностям не переходим
#Дики-Фуллера вполне достаточно,но можно еще проделать КПСС-тест и ПП-тест,они также покажут наличие стационарности
kpss.test(Sample1$Ряд3, null='Level') #противоречит остальным и говорит, что ряд не стационарный, ну и фиг с ним, можно забить
## Warning in kpss.test(Sample1$Ряд3, null = "Level"): p-value greater than
## printed p-value
##
## KPSS Test for Level Stationarity
##
## data: Sample1$Ряд3
## KPSS Level = 0.14616, Truncation lag parameter = 7, p-value = 0.1
pp.test(Sample1$Ряд3) #p-value=0.01 что меньше(или равно) 0.01 => Null отвергаем, ряд стационарный
## Warning in pp.test(Sample1$Ряд3): p-value smaller than printed p-value
##
## Phillips-Perron Unit Root Test
##
## data: Sample1$Ряд3
## Dickey-Fuller Z(alpha) = -1344.8, Truncation lag parameter = 7, p-value
## = 0.01
## alternative hypothesis: stationary
#Если ряд оказался нестационарным, то мы берем first difference, если и он окажется нестационарным, то берем second difference
#Делается это так: diff1_data<-diff(data, differences=1) и diff2_data<-diff(data, differences=2) соответственно
c <- Sample1$Ряд3
#Либо проверяем так:
aTSA::adf.test(c)
## Augmented Dickey-Fuller Test
## alternative: stationary
##
## Type 1: no drift no trend
## lag ADF p.value
## [1,] 0 -43.12 0.01
## [2,] 1 -21.96 0.01
## [3,] 2 -15.51 0.01
## [4,] 3 -12.69 0.01
## [5,] 4 -10.74 0.01
## [6,] 5 -9.30 0.01
## [7,] 6 -8.35 0.01
## Type 2: with drift no trend
## lag ADF p.value
## [1,] 0 -47.8 0.01
## [2,] 1 -25.9 0.01
## [3,] 2 -19.3 0.01
## [4,] 3 -16.6 0.01
## [5,] 4 -14.7 0.01
## [6,] 5 -13.3 0.01
## [7,] 6 -12.3 0.01
## Type 3: with drift and trend
## lag ADF p.value
## [1,] 0 -47.8 0.01
## [2,] 1 -25.9 0.01
## [3,] 2 -19.3 0.01
## [4,] 3 -16.6 0.01
## [5,] 4 -14.7 0.01
## [6,] 5 -13.3 0.01
## [7,] 6 -12.4 0.01
## ----
## Note: in fact, p.value = 0.01 means p.value <= 0.01
aTSA::kpss.test(c)
## KPSS Unit Root Test
## alternative: nonstationary
##
## Type 1: no drift no trend
## lag stat p.value
## 7 24.4 0.01
## -----
## Type 2: with drift no trend
## lag stat p.value
## 7 0.177 0.1
## -----
## Type 1: with drift and trend
## lag stat p.value
## 7 0.0266 0.1
## -----------
## Note: p.value = 0.01 means p.value <= 0.01
## : p.value = 0.10 means p.value >= 0.10
aTSA::pp.test(c)
## Phillips-Perron Unit Root Test
## alternative: stationary
##
## Type 1: no drift no trend
## lag Z_rho p.value
## 7 -1590 0.01
## -----
## Type 2: with drift no trend
## lag Z_rho p.value
## 7 -1347 0.01
## -----
## Type 3: with drift and trend
## lag Z_rho p.value
## 7 -1345 0.01
## ---------------
## Note: p-value = 0.01 means p.value <= 0.01
#По результатам тестов, ряд стационарный
#3.Структурные разрывы=========================================================================================
#Используем Sup-F тест в пакете strucchange
#Построим переменную с лагом 1
data1 <-c(0,Sample1$Ряд3[1:length(Sample1$Ряд3)-1])
#Считаем F-статистику
stat <- Fstats(Sample1$Ряд3~data1,from=0.1)
#Аlpha-уровень значимости для Чоу-теста
plot(stat,alpha=0.01)
#Определим точку разрыва
breakpoints(stat)
##
## Optimal 2-segment partition:
##
## Call:
## breakpoints.Fstats(obj = stat)
##
## Breakpoints at observation number:
## 100
##
## Corresponding to breakdates:
## 0.1035565
sctest(stat,type="supF")
##
## supF test
##
## data: stat
## sup.F = 9.8282, p-value = 0.1242
#Построим функцию Байесовского информационного критерия-БИК
stata <- breakpoints(Sample1$Ряд3~data1)
summary(stata)
##
## Optimal (m+1)-segment partition:
##
## Call:
## breakpoints.formula(formula = Sample1$Ряд3 ~ data1)
##
## Breakpoints at observation number:
##
## m = 1 365
## m = 2 365 638
## m = 3 365 638 809
## m = 4 218 365 638 809
## m = 5 174 347 504 647 809
##
## Corresponding to breakdates:
##
## m = 1 0.381799163179916
## m = 2 0.381799163179916 0.667364016736402
## m = 3 0.381799163179916 0.667364016736402
## m = 4 0.228033472803347 0.381799163179916 0.667364016736402
## m = 5 0.182008368200837 0.362970711297071 0.527196652719665 0.676778242677824
##
## m = 1
## m = 2
## m = 3 0.846234309623431
## m = 4 0.846234309623431
## m = 5 0.846234309623431
##
## Fit:
##
## m 0 1 2 3 4 5
## RSS 7159 7109 7071 7046 7025 7025
## BIC 4658 4672 4688 4705 4723 4743
plot(stata)
## compute breakdates corresponding to the
## breakpoints of minimum BIC segmentation
breakdates(stata) #Вывел NA - точек разрыва нет
## [1] NA
#4.Разбиваем ряд и добиваемся стационарности при необходимости
#Стационарным ряд у нас был изначально => к первым разностям мы не переходим
#Точек разрыва не имеет => ряд не разбиваем
#5.Строим ARMA-модель с учетом структурных разрывов(если они есть)
eacf(Sample1$Ряд3)
## AR/MA
## 0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 x x o o o o o o o o o o o o
## 1 o o o o o o o o o o o o o o
## 2 o o o o o o o o o o o o o o
## 3 x x o o o o o o o o o o o o
## 4 x o x o o o o o o o o o o o
## 5 x o o x o o o o o o o o o o
## 6 x x x x o o o o o o o o o o
## 7 x x x x x x o o o o o o o o
#Из матрицы EACF делаем вывод о том, что оптимальна модель arma(1,0)
# Оцениваем параметры оптимальной “короткой” модели (на основе логарифма индекса)
m=Arima(c, order=c(1,0,0))
coeftest(m)
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -0.411905 0.029538 -13.945 < 2.2e-16 ***
## intercept 0.779931 0.062674 12.444 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#6.Тестируем остатки модели на автокорреляцию. Корректируем модель, чтобы добиться отсутствия автокорреляции.
#Проверяем остатки на автокоррелированность (ACF,Ljung-Box test)
plot(residuals(m))
Acf(residuals(m))
#По графику видно, что автокорреляции остатков нет
Box.test(residuals(m), lag=6, fitdf=1, type = "Ljung-Box")
##
## Box-Ljung test
##
## data: residuals(m)
## X-squared = 1.1306, df = 5, p-value = 0.9514
#Null об отсутсвии в тесте Льюинга-Бокса автокорреляции не отвергаем, поскольку p-value = 0.96
#7.Тестируем наличие ARCH-эффекта в остатках.Строим ARMA + GARCH, если ARCH-эффект есть
Acf(residuals(m)^2)
Box.test(residuals(m)^2, lag=6, fitdf=1, type = "Ljung-Box")
##
## Box-Ljung test
##
## data: residuals(m)^2
## X-squared = 121.9, df = 5, p-value < 2.2e-16
#По графику видно, что автокорреляция квадратов остатков имеется, маленький p-value также отвергает гипотезу об отсутсвии автокорреляции в тесте Льюинга-Бокса
Pacf(residuals(m)^2)
eacf(residuals(m)^2)
## AR/MA
## 0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 x x x o o o o x o o x o o o
## 1 x x x x o o o o o o x o o o
## 2 x x x o o o o o o o x o o o
## 3 x x x o o o o o o o x o o x
## 4 x x x o x o o o o o x o x o
## 5 x x x x o o o o o o x o x o
## 6 x x x x x o o o o o x o x o
## 7 x x x x x o o o o o x o x o
#Строим стандартную гарч-модель (1,1)
spec = ugarchspec(variance.model = list(model = 'sGARCH',garchOrder = c(1, 1)),
mean.model = list(armaOrder = c(1,0), include.mean = TRUE), distribution.model = "norm")
m = ugarchfit(spec, c)
m
##
## *---------------------------------*
## * GARCH Model Fit *
## *---------------------------------*
##
## Conditional Variance Dynamics
## -----------------------------------
## GARCH Model : sGARCH(1,1)
## Mean Model : ARFIMA(1,0,0)
## Distribution : norm
##
## Optimal Parameters
## ------------------------------------
## Estimate Std. Error t value Pr(>|t|)
## mu 0.75063 0.056691 13.2408 0.000000
## ar1 -0.37880 0.032410 -11.6878 0.000000
## omega 0.88131 0.311291 2.8312 0.004638
## alpha1 0.19843 0.037637 5.2722 0.000000
## beta1 0.68903 0.063983 10.7689 0.000000
##
## Robust Standard Errors:
## Estimate Std. Error t value Pr(>|t|)
## mu 0.75063 0.061467 12.2119 0.000000
## ar1 -0.37880 0.037365 -10.1380 0.000000
## omega 0.88131 0.348238 2.5308 0.011381
## alpha1 0.19843 0.043283 4.5844 0.000005
## beta1 0.68903 0.072073 9.5602 0.000000
##
## LogLikelihood : -2272.493
##
## Information Criteria
## ------------------------------------
##
## Akaike 4.7646
## Bayes 4.7901
## Shibata 4.7646
## Hannan-Quinn 4.7743
##
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 0.005716 0.9397
## Lag[2*(p+q)+(p+q)-1][2] 0.264348 0.9965
## Lag[4*(p+q)+(p+q)-1][5] 1.115716 0.9281
## d.o.f=1
## H0 : No serial correlation
##
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 1.245 0.26453
## Lag[2*(p+q)+(p+q)-1][5] 4.474 0.20055
## Lag[4*(p+q)+(p+q)-1][9] 8.596 0.09819
## d.o.f=2
##
## Weighted ARCH LM Tests
## ------------------------------------
## Statistic Shape Scale P-Value
## ARCH Lag[3] 2.164 0.500 2.000 0.14130
## ARCH Lag[5] 6.728 1.440 1.667 0.04042
## ARCH Lag[7] 8.521 2.315 1.543 0.04013
##
## Nyblom stability test
## ------------------------------------
## Joint Statistic: 0.4717
## Individual Statistics:
## mu 0.14063
## ar1 0.08844
## omega 0.08091
## alpha1 0.10691
## beta1 0.06747
##
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic: 1.28 1.47 1.88
## Individual Statistic: 0.35 0.47 0.75
##
## Sign Bias Test
## ------------------------------------
## t-value prob sig
## Sign Bias 0.3864 0.6993
## Negative Sign Bias 1.0143 0.3107
## Positive Sign Bias 0.1122 0.9107
## Joint Effect 1.0791 0.7821
##
##
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
## group statistic p-value(g-1)
## 1 20 28.69 0.07108
## 2 30 31.24 0.35427
## 3 40 48.35 0.14490
## 4 50 55.92 0.23096
##
##
## Elapsed time : 0.09973311
Acf (residuals(m, standardize="TRUE"))
Box.test(residuals(m,standardize="TRUE"), lag=6, fitdf=1, type = "Ljung-Box")
##
## Box-Ljung test
##
## data: residuals(m, standardize = "TRUE")
## X-squared = 2.1256, df = 5, p-value = 0.8315
Acf (residuals(m, standardize="TRUE")^2)
Box.test(residuals(m,standardize="TRUE"), lag=6, fitdf=1, type = "Ljung-Box")
##
## Box-Ljung test
##
## data: residuals(m, standardize = "TRUE")
## X-squared = 2.1256, df = 5, p-value = 0.8315
plot(residuals(m,standardize="TRUE"))
#8.Запишем формулу итоговой модели
#Выписать estimate из выводов функций coeftest(m) и optimal parameters функции m = ugarchfit(spec, c) m
#-----------------------------------
#GARCH Model : sGARCH(1,1)
#Mean Model : ARFIMA(1,0,0)
#Distribution : norm
#z test of coefficients:
#------------------------------------
# Estimate Std. Error z value Pr(>|z|)
#ar1 -0.411905 0.029538 -13.945 < 2.2e-16 ***
#intercept 0.779931 0.062674 12.444 < 2.2e-16 ***
#Optimal Parameters
#------------------------------------
# Estimate Std. Error t value Pr(>|t|)
#mu 0.75063 0.056691 13.2408 0.000000
#ar1 -0.37880 0.032410 -11.6878 0.000000
#omega 0.88131 0.311291 2.8312 0.004638
#alpha1 0.19843 0.037637 5.2722 0.000000
#beta1 0.68903 0.063983 10.7689 0.000000