library(quantmod)
## Загрузка требуемого пакета: xts
## Загрузка требуемого пакета: zoo
## 
## Присоединяю пакет: 'zoo'
## Следующие объекты скрыты от 'package:base':
## 
##     as.Date, as.Date.numeric
## Загрузка требуемого пакета: TTR
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library(forecast)
library(tseries)
library(rugarch)
## Загрузка требуемого пакета: parallel
## 
## Присоединяю пакет: 'rugarch'
## Следующий объект скрыт от 'package:stats':
## 
##     sigma
library(ggplot2)
library(stargazer)
## 
## Please cite as:
##  Hlavac, Marek (2022). stargazer: Well-Formatted Regression and Summary Statistics Tables.
##  R package version 5.2.3. https://CRAN.R-project.org/package=stargazer
library(lmtest)
library(tseries)
library(urca)
library(FinTS)
## 
## Присоединяю пакет: 'FinTS'
## Следующий объект скрыт от 'package:forecast':
## 
##     Acf
SIRI <- getSymbols("SIRI", from = "2010-01-01", to = "2022-12-31", auto.assign = FALSE)
Price <- SIRI$SIRI.Close
chartSeries(SIRI)

adf.test(Price)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  Price
## Dickey-Fuller = -3.0082, Lag order = 14, p-value = 0.1514
## alternative hypothesis: stationary
#Стационарны ли первые разности? 
adf.test(na.omit(diff(Price))) #ряд стационарен
## Warning in adf.test(na.omit(diff(Price))): p-value smaller than printed p-value
## 
##  Augmented Dickey-Fuller Test
## 
## data:  na.omit(diff(Price))
## Dickey-Fuller = -15.827, Lag order = 14, p-value = 0.01
## alternative hypothesis: stationary
summary(ur.df(na.omit(diff(Price))))
## 
## ############################################### 
## # 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 
## -0.78135 -0.03273  0.00258  0.03870  0.49704 
## 
## Coefficients:
##            Estimate Std. Error t value Pr(>|t|)    
## z.lag.1    -1.04272    0.02563 -40.678   <2e-16 ***
## z.diff.lag -0.02906    0.01749  -1.662   0.0967 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.0795 on 3267 degrees of freedom
## Multiple R-squared:  0.5374, Adjusted R-squared:  0.5371 
## F-statistic:  1897 on 2 and 3267 DF,  p-value: < 2.2e-16
## 
## 
## Value of test-statistic is: -40.6782 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau1 -2.58 -1.95 -1.62
dim(SIRI$SIRI.Close %>% diff())
## [1] 3272    1
SIRI$SIRI.Close %>% diff() %>% ggtsdisplay(main="", theme = theme_classic()) 
## Warning: Removed 1 rows containing missing values (`geom_point()`).

SIRI$SIRI.Close %>% diff() %>% ggtsdisplay(main="", theme = theme_classic(), plot.type="scatter") 
## Warning: Removed 1 rows containing missing values (`geom_point()`).
## Warning: Removed 1 rows containing missing values (`geom_point()`).

#значения ACF и PACF
acf(na.omit(SIRI$SIRI.Close %>% diff()), plot = FALSE)
## 
## Autocorrelations of series 'na.omit(SIRI$SIRI.Close %>% diff())', by lag
## 
##      0      1      2      3      4      5      6      7      8      9     10 
##  1.000 -0.074  0.034 -0.014 -0.005  0.029 -0.058  0.031 -0.075  0.004  0.009 
##     11     12     13     14     15     16     17     18     19     20     21 
## -0.012  0.016 -0.021  0.007 -0.025 -0.018 -0.002 -0.011  0.001 -0.010  0.018 
##     22     23     24     25     26     27     28     29     30     31     32 
## -0.034  0.017  0.011 -0.002 -0.012 -0.002 -0.015 -0.022 -0.001 -0.048 -0.004 
##     33     34     35 
##  0.006 -0.006  0.000
pacf(na.omit(SIRI$SIRI.Close %>% diff()), plot = FALSE)
## 
## Partial autocorrelations of series 'na.omit(SIRI$SIRI.Close %>% diff())', by lag
## 
##      1      2      3      4      5      6      7      8      9     10     11 
## -0.074  0.029 -0.009 -0.007  0.029 -0.054  0.021 -0.068 -0.009  0.012 -0.009 
##     12     13     14     15     16     17     18     19     20     21     22 
##  0.009 -0.013 -0.004 -0.021 -0.026 -0.007 -0.007 -0.005 -0.006  0.013 -0.034 
##     23     24     25     26     27     28     29     30     31     32     33 
##  0.009  0.011 -0.001 -0.016  0.000 -0.021 -0.022 -0.008 -0.049 -0.010  0.006 
##     34     35 
## -0.008 -0.005
#построение моделей
mod1 <- auto.arima(Price) #модель - ARIMA(0,1,2)
mod2 <- arima(Price, order = c(0,1,1))
mod3 <- arima(Price, order = c(0,1,2)) #для сохранения через stargazer
mod4 <- arima(Price, order = c(0,1,3))

mod1
## Series: Price 
## ARIMA(0,1,2) 
## 
## Coefficients:
##           ma1     ma2
##       -0.0710  0.0324
## s.e.   0.0175  0.0174
## 
## sigma^2 = 0.006318:  log likelihood = 3642.49
## AIC=-7278.98   AICc=-7278.97   BIC=-7260.7
mod2
## 
## Call:
## arima(x = Price, order = c(0, 1, 1))
## 
## Coefficients:
##           ma1
##       -0.0697
## s.e.   0.0169
## 
## sigma^2 estimated as 0.00632:  log likelihood = 3640.75,  aic = -7277.51
mod3
## 
## Call:
## arima(x = Price, order = c(0, 1, 2))
## 
## Coefficients:
##           ma1     ma2
##       -0.0710  0.0324
## s.e.   0.0175  0.0174
## 
## sigma^2 estimated as 0.006314:  log likelihood = 3642.49,  aic = -7278.98
mod4
## 
## Call:
## arima(x = Price, order = c(0, 1, 3))
## 
## Coefficients:
##           ma1     ma2      ma3
##       -0.0714  0.0341  -0.0164
## s.e.   0.0175  0.0175   0.0184
## 
## sigma^2 estimated as 0.006312:  log likelihood = 3642.89,  aic = -7277.77
AIC(mod1)
## [1] -7278.978
AIC(mod2)
## [1] -7277.509
#первые 2000 наблюдений меньше волатильности, чем на последних 1000
Price1 <- na.omit(SIRI$SIRI.Close %>% diff())
var.test(SIRI$SIRI.Close[c(1:2000),], SIRI$SIRI.Close[c(2000:3000),]) #да это так
## 
##  F test to compare two variances
## 
## data:  SIRI$SIRI.Close[c(1:2000), ] and SIRI$SIRI.Close[c(2000:3000), ]
## F = 5.9722, num df = 1999, denom df = 1000, p-value < 2.2e-16
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  5.358728 6.642889
## sample estimates:
##            SIRI.Close
## SIRI.Close   5.972171
## attr(,"names")
## [1] "ratio of variances"
var.test(Price1[c(1:2000),], Price1[c(2000:3000),])
## 
##  F test to compare two variances
## 
## data:  Price1[c(1:2000), ] and Price1[c(2000:3000), ]
## F = 0.1937, num df = 1999, denom df = 1000, p-value < 2.2e-16
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.1738000 0.2154493
## sample estimates:
##            SIRI.Close
## SIRI.Close  0.1936959
## attr(,"names")
## [1] "ratio of variances"
#проверка на автокорреляцию
checkresiduals(mod3)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(0,1,2)
## Q* = 30.572, df = 8, p-value = 0.0001674
## 
## Model df: 2.   Total lags used: 10
#видна автокорреляция, кластеры волатильности

#есть ли ARCH структара?
ArchTest(resid(mod3))
## 
##  ARCH LM-test; Null hypothesis: no ARCH effects
## 
## data:  resid(mod3)
## Chi-squared = 564.34, df = 12, p-value < 2.2e-16
#да

#нормальны ли остатки?
jarque.bera.test(resid(mod3))
## 
##  Jarque Bera Test
## 
## data:  resid(mod3)
## X-squared = 19991, df = 2, p-value < 2.2e-16
#нет

#попробуем добавить гетероскедтичность
e <- resid(mod3)
acf(e^2)

pacf(e^2)

#похоже на GARCH структуру
e^2 %>% ggtsdisplay(main="", theme = theme_classic()) 

#GARCH модель
garch_spec1 <- ugarchspec(variance.model=list(model="sGARCH", garchOrder=c(1,1)), mean.model=list(armaOrder=c(0,2)))
fit_garch1 <- ugarchfit(spec = garch_spec1, data = na.omit(diff(Price)))
fit_garch1
## 
## *---------------------------------*
## *          GARCH Model Fit        *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## -----------------------------------
## GARCH Model  : sGARCH(1,1)
## Mean Model   : ARFIMA(0,0,2)
## Distribution : norm 
## 
## Optimal Parameters
## ------------------------------------
##         Estimate  Std. Error  t value Pr(>|t|)
## mu      0.002258    0.000948   2.3815 0.017244
## ma1    -0.033328    0.018349  -1.8164 0.069311
## ma2    -0.018796    0.018573  -1.0120 0.311515
## omega   0.000015    0.000005   2.9999 0.002700
## alpha1  0.037831    0.002614  14.4715 0.000000
## beta1   0.961169    0.001844 521.2528 0.000000
## 
## Robust Standard Errors:
##         Estimate  Std. Error   t value Pr(>|t|)
## mu      0.002258    0.000918   2.46093 0.013858
## ma1    -0.033328    0.019796  -1.68356 0.092267
## ma2    -0.018796    0.020683  -0.90878 0.363464
## omega   0.000015    0.000011   1.35050 0.176857
## alpha1  0.037831    0.005945   6.36334 0.000000
## beta1   0.961169    0.002621 366.74926 0.000000
## 
## LogLikelihood : 4213.048 
## 
## Information Criteria
## ------------------------------------
##                     
## Akaike       -2.5723
## Bayes        -2.5612
## Shibata      -2.5723
## Hannan-Quinn -2.5683
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                  0.0005998  0.9805
## Lag[2*(p+q)+(p+q)-1][5] 0.2048564  1.0000
## Lag[4*(p+q)+(p+q)-1][9] 1.9190880  0.9877
## d.o.f=2
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                      3.169 0.07506
## Lag[2*(p+q)+(p+q)-1][5]     6.047 0.08770
## Lag[4*(p+q)+(p+q)-1][9]     6.826 0.21419
## d.o.f=2
## 
## Weighted ARCH LM Tests
## ------------------------------------
##             Statistic Shape Scale P-Value
## ARCH Lag[3]    0.1833 0.500 2.000  0.6686
## ARCH Lag[5]    0.3019 1.440 1.667  0.9398
## ARCH Lag[7]    0.7671 2.315 1.543  0.9482
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  0.9908
## Individual Statistics:              
## mu     0.11258
## ma1    0.09702
## ma2    0.03604
## omega  0.21126
## alpha1 0.27390
## beta1  0.25369
## 
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic:          1.49 1.68 2.12
## Individual Statistic:     0.35 0.47 0.75
## 
## Sign Bias Test
## ------------------------------------
##                    t-value      prob sig
## Sign Bias           2.1107 0.0348732  **
## Negative Sign Bias  0.5489 0.5830977    
## Positive Sign Bias  3.3564 0.0007986 ***
## Joint Effect       12.2711 0.0065099 ***
## 
## 
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
##   group statistic p-value(g-1)
## 1    20     198.0    8.409e-32
## 2    30     321.0    5.564e-51
## 3    40     307.8    1.725e-43
## 4    50     315.9    1.090e-40
## 
## 
## Elapsed time : 0.328692
garch_spec2 <- ugarchspec(variance.model=list(model="sGARCH", garchOrder=c(1,2)), mean.model=list(armaOrder=c(0,2)))
fit_garch2 <- ugarchfit(spec = garch_spec2, data = na.omit(diff(Price)))
fit_garch2
## 
## *---------------------------------*
## *          GARCH Model Fit        *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## -----------------------------------
## GARCH Model  : sGARCH(1,2)
## Mean Model   : ARFIMA(0,0,2)
## Distribution : norm 
## 
## Optimal Parameters
## ------------------------------------
##         Estimate  Std. Error  t value Pr(>|t|)
## mu      0.002230    0.000944   2.3614 0.018205
## ma1    -0.033175    0.018692  -1.7748 0.075937
## ma2    -0.019223    0.018372  -1.0463 0.295436
## omega   0.000019    0.000008   2.4538 0.014137
## alpha1  0.053922    0.007111   7.5831 0.000000
## beta1   0.533068    0.015769  33.8042 0.000000
## beta2   0.412010    0.015189  27.1263 0.000000
## 
## Robust Standard Errors:
##         Estimate  Std. Error  t value Pr(>|t|)
## mu      0.002230    0.000917  2.43159 0.015033
## ma1    -0.033175    0.019659 -1.68749 0.091510
## ma2    -0.019223    0.019915 -0.96522 0.334435
## omega   0.000019    0.000020  0.94823 0.343011
## alpha1  0.053922    0.026360  2.04557 0.040798
## beta1   0.533068    0.014061 37.91068 0.000000
## beta2   0.412010    0.013269 31.05067 0.000000
## 
## LogLikelihood : 4216.169 
## 
## Information Criteria
## ------------------------------------
##                     
## Akaike       -2.5736
## Bayes        -2.5606
## Shibata      -2.5736
## Hannan-Quinn -2.5690
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                  0.0008563  0.9767
## Lag[2*(p+q)+(p+q)-1][5] 0.2239654  1.0000
## Lag[4*(p+q)+(p+q)-1][9] 1.9154667  0.9878
## d.o.f=2
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                          statistic p-value
## Lag[1]                       1.700  0.1923
## Lag[2*(p+q)+(p+q)-1][8]      5.740  0.2622
## Lag[4*(p+q)+(p+q)-1][14]     6.875  0.5221
## d.o.f=3
## 
## Weighted ARCH LM Tests
## ------------------------------------
##             Statistic Shape Scale P-Value
## ARCH Lag[4]    0.1487 0.500 2.000  0.6997
## ARCH Lag[6]    0.2904 1.461 1.711  0.9476
## ARCH Lag[8]    0.7675 2.368 1.583  0.9553
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  1.0935
## Individual Statistics:              
## mu     0.11142
## ma1    0.08874
## ma2    0.03336
## omega  0.43602
## alpha1 0.25858
## beta1  0.24545
## beta2  0.24304
## 
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic:          1.69 1.9 2.35
## Individual Statistic:     0.35 0.47 0.75
## 
## Sign Bias Test
## ------------------------------------
##                    t-value     prob sig
## Sign Bias           2.0879 0.036887  **
## Negative Sign Bias  0.2337 0.815249    
## Positive Sign Bias  3.0325 0.002445 ***
## Joint Effect        9.9740 0.018788  **
## 
## 
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
##   group statistic p-value(g-1)
## 1    20     203.1    8.319e-33
## 2    30     320.2    7.987e-51
## 3    40     304.1    8.667e-43
## 4    50     322.6    6.172e-42
## 
## 
## Elapsed time : 0.333338
garch_spec3 <- ugarchspec(variance.model=list(model="sGARCH", garchOrder=c(2,2)), mean.model=list(armaOrder=c(0,2)))
fit_garch3 <- ugarchfit(spec = garch_spec3, data = na.omit(diff(Price)))
fit_garch3
## 
## *---------------------------------*
## *          GARCH Model Fit        *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## -----------------------------------
## GARCH Model  : sGARCH(2,2)
## Mean Model   : ARFIMA(0,0,2)
## Distribution : norm 
## 
## Optimal Parameters
## ------------------------------------
##         Estimate  Std. Error   t value Pr(>|t|)
## mu      0.002230    0.000945  2.361068 0.018222
## ma1    -0.033158    0.018692 -1.773860 0.076086
## ma2    -0.019217    0.018373 -1.045949 0.295585
## omega   0.000019    0.000008  2.380103 0.017308
## alpha1  0.053915    0.013367  4.033564 0.000055
## alpha2  0.000000    0.013625  0.000001 0.999999
## beta1   0.533382    0.015707 33.958187 0.000000
## beta2   0.411703    0.014776 27.863772 0.000000
## 
## Robust Standard Errors:
##         Estimate  Std. Error   t value Pr(>|t|)
## mu      0.002230    0.000913  2.443937 0.014528
## ma1    -0.033158    0.019670 -1.685737 0.091846
## ma2    -0.019217    0.019920 -0.964707 0.334692
## omega   0.000019    0.000022  0.891565 0.372626
## alpha1  0.053915    0.023910  2.254918 0.024138
## alpha2  0.000000    0.021766  0.000001 0.999999
## beta1   0.533382    0.013326 40.026957 0.000000
## beta2   0.411703    0.016010 25.715706 0.000000
## 
## LogLikelihood : 4216.169 
## 
## Information Criteria
## ------------------------------------
##                     
## Akaike       -2.5730
## Bayes        -2.5581
## Shibata      -2.5730
## Hannan-Quinn -2.5677
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                  0.0008092  0.9773
## Lag[2*(p+q)+(p+q)-1][5] 0.2236244  1.0000
## Lag[4*(p+q)+(p+q)-1][9] 1.9151226  0.9878
## d.o.f=2
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                          statistic p-value
## Lag[1]                       1.700  0.1923
## Lag[2*(p+q)+(p+q)-1][11]     6.187  0.4110
## Lag[4*(p+q)+(p+q)-1][19]     7.838  0.6938
## d.o.f=4
## 
## Weighted ARCH LM Tests
## ------------------------------------
##             Statistic Shape Scale P-Value
## ARCH Lag[5]   0.05333 0.500 2.000  0.8174
## ARCH Lag[7]   0.63293 1.473 1.746  0.8607
## ARCH Lag[9]   0.85835 2.402 1.619  0.9495
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  3.8964
## Individual Statistics:              
## mu     0.11135
## ma1    0.08858
## ma2    0.03335
## omega  0.43500
## alpha1 0.25827
## alpha2 0.30143
## beta1  0.24514
## beta2  0.24273
## 
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic:          1.89 2.11 2.59
## Individual Statistic:     0.35 0.47 0.75
## 
## Sign Bias Test
## ------------------------------------
##                    t-value     prob sig
## Sign Bias           2.0878 0.036890  **
## Negative Sign Bias  0.2336 0.815330    
## Positive Sign Bias  3.0323 0.002446 ***
## Joint Effect        9.9729 0.018798  **
## 
## 
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
##   group statistic p-value(g-1)
## 1    20     202.9    8.848e-33
## 2    30     320.0    8.614e-51
## 3    40     303.8    1.018e-42
## 4    50     322.0    8.016e-42
## 
## 
## Elapsed time : 0.4039321
infocriteria(fit_garch1)
##                       
## Akaike       -2.572331
## Bayes        -2.561155
## Shibata      -2.572338
## Hannan-Quinn -2.568329
infocriteria(fit_garch2)
##                       
## Akaike       -2.573628
## Bayes        -2.560589
## Shibata      -2.573637
## Hannan-Quinn -2.568959
infocriteria(fit_garch3)
##                       
## Akaike       -2.573017
## Bayes        -2.558115
## Shibata      -2.573029
## Hannan-Quinn -2.567680
infocriteria(fit_garch1) > infocriteria(fit_garch2)
##                   
## Akaike        TRUE
## Bayes        FALSE
## Shibata       TRUE
## Hannan-Quinn  TRUE
infocriteria(fit_garch1) > infocriteria(fit_garch3)
##                   
## Akaike        TRUE
## Bayes        FALSE
## Shibata       TRUE
## Hannan-Quinn FALSE
infocriteria(fit_garch2) < infocriteria(fit_garch3)
##                  
## Akaike       TRUE
## Bayes        TRUE
## Shibata      TRUE
## Hannan-Quinn TRUE
#AIC меньший в модели fit_garch2

#визуальная диагностика
#прогнозы: 2SD
plot(fit_garch2, which = 1)

#VAR
plot(fit_garch2, which = 2)
## 
## please wait...calculating quantiles...

#условная дисперсия
plot(fit_garch2, which = 3)

#распределение остатков
plot(fit_garch2, which = 8)