Promoting the use of R in the NHS

Blog Article

This post was originally published on this site

(This article was first published on R Programming – DataScience+, and kindly contributed to R-bloggers)

    Categories

    1. Advanced Modeling

    Tags

    1. Data Visualisation
    2. Linear Regression
    3. R Programming

    This is the final part of the 4-series posts. In this fourth post, I am going to build an ARMA-GARCH model for Dow Jones Industrial Average (DJIA) daily trade volume log ratio. You can read the other three parts in the following links: part 1, part2, and part 3.

    Packages

    The packages being used in this post series are herein listed.

    suppressPackageStartupMessages(library(lubridate))
    suppressPackageStartupMessages(library(fBasics))
    suppressPackageStartupMessages(library(lmtest))
    suppressPackageStartupMessages(library(urca))
    suppressPackageStartupMessages(library(ggplot2))
    suppressPackageStartupMessages(library(quantmod))
    suppressPackageStartupMessages(library(PerformanceAnalytics))
    suppressPackageStartupMessages(library(rugarch))
    suppressPackageStartupMessages(library(FinTS))
    suppressPackageStartupMessages(library(forecast))
    suppressPackageStartupMessages(library(strucchange))
    suppressPackageStartupMessages(library(TSA))
    

    Getting Data

    We upload the environment status as saved at the end of part 3.

    load(file='DowEnvironment.RData')
    

    Structural Breaks in Daily Trade Volume

    In part 2, we already highlighted some level shifts occurring within the daily trade volume. Here is the plot we already commented in part 2.

    plot(dj_vol)
    

    Now we want to identify and date such level shifts. First, we verify that a linear regression with constant mean is statistically significative.

    lm_dj_vol 
    ## 
    ## Call:
    ## lm(formula = dj_vol ~ 1, data = dj_vol)
    ## 
    ## Residuals:
    ##        Min         1Q     Median         3Q        Max 
    ## -193537268  -91369768  -25727268   65320232  698562732 
    ## 
    ## Coefficients:
    ##              Estimate Std. Error t value Pr(>|t|)    
    ## (Intercept) 201947268    2060772      98   

    Such linear regression model is statistically significative, hence it makes sense to proceed on with the breakpoints structural identification with respect a constant value.

    bp_dj_vol 
    ## 
    ##   Optimal (m+1)-segment partition: 
    ## 
    ## Call:
    ## breakpoints.formula(formula = dj_vol ~ 1, h = 0.1)
    ## 
    ## Breakpoints at observation number:
    ##                                             
    ## m = 1                                   2499
    ## m = 2           896                     2499
    ## m = 3       626     1254                2499
    ## m = 4   342 644     1254                2499
    ## m = 5   342 644     1219 1649           2499
    ## m = 6   320 622 924 1251 1649           2499
    ## m = 7   320 622 924 1251 1692      2172 2499
    ## m = 8   320 622 924 1251 1561 1863 2172 2499
    ## 
    ## Corresponding to breakdates:
    ##                                                              
    ## m = 1                                                        
    ## m = 2                                       0.296688741721854
    ## m = 3                     0.207284768211921                  
    ## m = 4   0.113245033112583 0.213245033112583                  
    ## m = 5   0.113245033112583 0.213245033112583                  
    ## m = 6   0.105960264900662 0.205960264900662 0.305960264900662
    ## m = 7   0.105960264900662 0.205960264900662 0.305960264900662
    ## m = 8   0.105960264900662 0.205960264900662 0.305960264900662
    ##                                                              
    ## m = 1                                                        
    ## m = 2                                                        
    ## m = 3   0.41523178807947                                     
    ## m = 4   0.41523178807947                                     
    ## m = 5   0.40364238410596  0.546026490066225                  
    ## m = 6   0.414238410596027 0.546026490066225                  
    ## m = 7   0.414238410596027 0.560264900662252                  
    ## m = 8   0.414238410596027 0.516887417218543 0.616887417218543
    ##                                            
    ## m = 1                     0.827483443708609
    ## m = 2                     0.827483443708609
    ## m = 3                     0.827483443708609
    ## m = 4                     0.827483443708609
    ## m = 5                     0.827483443708609
    ## m = 6                     0.827483443708609
    ## m = 7   0.719205298013245 0.827483443708609
    ## m = 8   0.719205298013245 0.827483443708609
    ## 
    ## Fit:
    ##                                                                          
    ## m   0         1         2         3         4         5         6        
    ## RSS 3.872e+19 2.772e+19 1.740e+19 1.547e+19 1.515e+19 1.490e+19 1.475e+19
    ## BIC 1.206e+05 1.196e+05 1.182e+05 1.179e+05 1.178e+05 1.178e+05 1.178e+05
    ##                        
    ## m   7         8        
    ## RSS 1.472e+19 1.478e+19
    ## BIC 1.178e+05 1.178e+05
    

    It gives this plot.

    plot(bp_dj_vol)
    

    The minimum BIC is reached at breaks = 6.

    Here is the plot of the DJIA daily trade volume with level shifts (red line).

    bp_fit 

    Here are the calendar dates when such level shifts occur.

    bd 
    ## [1] "2008-04-10" "2009-06-22" "2010-09-01" "2011-12-16" "2013-07-22"
    ## [6] "2016-12-02"
    

    See ref. [10] for further details about structural breaks.

    Daily Trade Volume Log Ratio Model

    As shown in part 2, the daily trade volume log ratio:

    [
    v_{t} := ln frac{V_{t}}{V_{t-1}}
    ]

    is already available within our data environment. We plot it.

    plot(dj_vol_log_ratio)
    

    Outliers Detection

    The Return.clean function within Performance Analytics package is able to clean return time series from outliers. Here below we compare the original time series with the outliers adjusted one.

    dj_vol_log_ratio_outliersadj 

    The prosecution of the analysis will be carried on with the original time series as a more conservative approach to volatility evaluation.

    Correlation plots

    Here below are the total and partial correlation plots.

    acf(dj_vol_log_ratio)
    

    pacf(dj_vol_log_ratio)
    

    Above plots may suggest some ARMA(p,q) model with p and q > 0. That will be verified within the prosecution of the present analysis.

    Unit root tests

    We run the Augmented Dickey-Fuller test as available within the urca package. The no trend and no drift test flavor is run.

    (urdftest_lag = floor(12* (nrow(dj_vol_log_ratio)/100)^0.25))
    ## [1] 28
    
    summary(ur.df(dj_vol_log_ratio, type = "none", lags = urdftest_lag, 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 
    ## -2.51742 -0.13544 -0.02485  0.12006  1.98641 
    ## 
    ## Coefficients:
    ##              Estimate Std. Error t value Pr(>|t|)    
    ## z.lag.1      -5.21571    0.25218 -20.682  

    Based on reported test statistics compared with critical values, we reject the null hypothesis of unit root presence. See ref. [6] for further details about the Augmented Dickey-Fuller test.

    ARMA model

    We now determine the ARMA structure of our time series in order to run the ARCH effects test on the resulting residuals. That is in agreement with what outlined in ref. [4] $4.3.

    We take advantage of auto.arima() function within forecast package (ref. [7]) to have an idea of what ARMA model to start with.

    (auto_model_1 
    ## Series: dj_vol_log_ratio 
    ## ARIMA(3,0,2) with zero mean 
    ## 
    ## Coefficients:
    ##           ar1     ar2     ar3      ma1      ma2
    ##       -0.5266  0.3769  0.1459  -0.0896  -0.7812
    ## s.e.   0.0970  0.0401  0.0203   0.0968   0.0872
    ## 
    ## sigma^2 estimated as 0.0659:  log likelihood=-176.48
    ## AIC=364.96   AICc=364.99   BIC=401.03
    
    coeftest(auto_model_1)
    ## 
    ## z test of coefficients:
    ## 
    ##      Estimate Std. Error z value  Pr(>|z|)    
    ## ar1 -0.526617   0.096988 -5.4297 5.644e-08 ***
    ## ar2  0.376906   0.040090  9.4014 

    The ma1 coefficient is not significative. Hence we try with the following ARMA(2,3) model.

    (arma_model_2 
    ## 
    ## Call:
    ## arima(x = dj_vol_log_ratio, order = c(2, 0, 3), include.mean = FALSE)
    ## 
    ## Coefficients:
    ##           ar1     ar2      ma1      ma2     ma3
    ##       -0.1802  0.6441  -0.4351  -0.8604  0.3596
    ## s.e.   0.0643  0.0454   0.0681   0.0493  0.0423
    ## 
    ## sigma^2 estimated as 0.0658:  log likelihood = -176.9,  aic = 363.79
    
    coeftest(arma_model_2)
    ## 
    ## z test of coefficients:
    ## 
    ##      Estimate Std. Error  z value  Pr(>|z|)    
    ## ar1 -0.180233   0.064315  -2.8023  0.005073 ** 
    ## ar2  0.644104   0.045449  14.1721 

    All coefficients are significative and AIC is lower than the one of the first model. We then try with ARMA(1,2).

    (arma_model_3 
    ## 
    ## Call:
    ## arima(x = dj_vol_log_ratio, order = c(1, 0, 2), include.mean = FALSE)
    ## 
    ## Coefficients:
    ##          ar1      ma1     ma2
    ##       0.6956  -1.3183  0.3550
    ## s.e.  0.0439   0.0518  0.0453
    ## 
    ## sigma^2 estimated as 0.06598:  log likelihood = -180.92,  aic = 367.84
    
    coeftest(arma_model_3)
    ## 
    ## z test of coefficients:
    ## 
    ##      Estimate Std. Error  z value  Pr(>|z|)    
    ## ar1  0.695565   0.043874  15.8537 

    This model has the highest AIC among the set and all coefficients are significative.

    We can also have a try with the eacf() function within the TSA package as further verification.

    eacf(dj_vol_log_ratio)
    ## AR/MA
    ##   0 1 2 3 4 5 6 7 8 9 10 11 12 13
    ## 0 x o o x x o o x o o x  o  o  o 
    ## 1 x x o x o o o x o o x  o  o  o 
    ## 2 x x x x o o o o o o x  o  o  o 
    ## 3 x x x x o o o o o o x  o  o  o 
    ## 4 x x x x x o o o o o x  o  o  o 
    ## 5 x x x x o o o o o o o  o  o  o 
    ## 6 x x x x x o x o o o o  o  o  o 
    ## 7 x x x x x o o o o o o  o  o  o
    

    The upper left triangle with “O” as a vertex seems to be located somehow within {(1,2),(2,2),(1,3),(2,3)}, which represents the set of potential (p,q) values according to each function output. To remark that we prefer to consider parsimonious models, that is why we do not go too far as AR and MA orders.

    We already verified ARMA models with (p,q) orders within the set {(3,2)(2,3)(1,2)}. Let us try {(2,2)(1,3)}

    (arma_model_4 
    ## 
    ## Call:
    ## arima(x = dj_vol_log_ratio, order = c(2, 0, 2), include.mean = FALSE)
    ## 
    ## Coefficients:
    ##          ar1      ar2      ma1     ma2
    ##       0.7174  -0.0096  -1.3395  0.3746
    ## s.e.  0.1374   0.0560   0.1361  0.1247
    ## 
    ## sigma^2 estimated as 0.06598:  log likelihood = -180.9,  aic = 369.8
    
    coeftest(arma_model_4)
    ## 
    ## z test of coefficients:
    ## 
    ##       Estimate Std. Error z value  Pr(>|z|)    
    ## ar1  0.7173631  0.1374135  5.2205 1.785e-07 ***
    ## ar2 -0.0096263  0.0560077 -0.1719  0.863536    
    ## ma1 -1.3394720  0.1361208 -9.8403 

    The ar2 coefficient is not significative.

    (arma_model_5 
    ## 
    ## Call:
    ## arima(x = dj_vol_log_ratio, order = c(1, 0, 3), include.mean = FALSE)
    ## 
    ## Coefficients:
    ##          ar1      ma1     ma2     ma3
    ##       0.7031  -1.3253  0.3563  0.0047
    ## s.e.  0.0657   0.0684  0.0458  0.0281
    ## 
    ## sigma^2 estimated as 0.06598:  log likelihood = -180.9,  aic = 369.8
    
    coeftest(arma_model_5)
    ## 
    ## z test of coefficients:
    ## 
    ##       Estimate Std. Error  z value  Pr(>|z|)    
    ## ar1  0.7030934  0.0656902  10.7032 

    The ma3 coefficient is not significative.

    As a conclusion, we keep both ARMA(1,2) and ARMA(2,3) as tentative mean models. We can now proceed on with ARCH effect tests.

    ARCH effect Test

    If ARCH effects are statistical significative for the residuals of our time series, a GARCH model is required.

    We test the candidate mean model ARMA(2,3).

    resid_dj_vol_log_ratio 
    ## 
    ##  ARCH LM-test; Null hypothesis: no ARCH effects
    ## 
    ## data:  resid_dj_vol_log_ratio - mean(resid_dj_vol_log_ratio)
    ## Chi-squared = 78.359, df = 12, p-value = 8.476e-12
    

    Based on reported p-value, we reject the null hypothesis of no ARCH effects.

    Let us have a look at the residual correlation plot too.

    par(mfrow=c(1,2))
    acf(resid_dj_vol_log_ratio)
    pacf(resid_dj_vol_log_ratio)
    

    We test the second candidate mean model, ARMA(1,2).

    resid_dj_vol_log_ratio 
    ## 
    ##  ARCH LM-test; Null hypothesis: no ARCH effects
    ## 
    ## data:  resid_dj_vol_log_ratio - mean(resid_dj_vol_log_ratio)
    ## Chi-squared = 74.768, df = 12, p-value = 4.065e-11
    

    Based on reported p-value, we reject the null hypothesis of no ARCH effects.

    Let us have a look at the residual correlation plots.

    par(mfrow=c(1,2))
    acf(resid_dj_vol_log_ratio)
    pacf(resid_dj_vol_log_ratio)
    

    To inspect asymmetries within the DJIA volume log ratio, summary statistics and density plot are shown.

    basicStats(dj_vol_log_ratio)
    ##              DJI.Volume
    ## nobs        3019.000000
    ## NAs            0.000000
    ## Minimum       -2.301514
    ## Maximum        2.441882
    ## 1. Quartile   -0.137674
    ## 3. Quartile    0.136788
    ## Mean          -0.000041
    ## Median        -0.004158
    ## Sum           -0.124733
    ## SE Mean        0.005530
    ## LCL Mean      -0.010885
    ## UCL Mean       0.010802
    ## Variance       0.092337
    ## Stdev          0.303869
    ## Skewness      -0.182683
    ## Kurtosis       9.463384
    

    Negative skewness is reported.

    plot(density(dj_vol_log_ratio))
    

    As a result, also for daily trade volume log ratio the eGARCH model will be proposed.

    We run two fits in order to compare the results with two candidate mean models ARMA(1,2) and ARMA(2,3)

    ARMA-GARCH: ARMA(1,2) + eGARCH(1,1)

    garchspec 
    ## 
    ## *---------------------------------*
    ## *          GARCH Model Fit        *
    ## *---------------------------------*
    ## 
    ## Conditional Variance Dynamics    
    ## -----------------------------------
    ## GARCH Model  : eGARCH(1,1)
    ## Mean Model   : ARFIMA(1,0,2)
    ## Distribution : sstd 
    ## 
    ## Optimal Parameters
    ## ------------------------------------
    ##         Estimate  Std. Error     t value Pr(>|t|)
    ## ar1      0.67731    0.014856     45.5918  0.0e+00
    ## ma1     -1.22817    0.000038 -31975.1819  0.0e+00
    ## ma2      0.27070    0.000445    608.3525  0.0e+00
    ## omega   -1.79325    0.207588     -8.6385  0.0e+00
    ## alpha1   0.14348    0.032569      4.4053  1.1e-05
    ## beta1    0.35819    0.073164      4.8957  1.0e-06
    ## gamma1   0.41914    0.042252      9.9199  0.0e+00
    ## skew     1.32266    0.031528     41.9518  0.0e+00
    ## shape    3.54346    0.221750     15.9795  0.0e+00
    ## 
    ## Robust Standard Errors:
    ##         Estimate  Std. Error     t value Pr(>|t|)
    ## ar1      0.67731    0.022072     30.6859  0.0e+00
    ## ma1     -1.22817    0.000067 -18466.0626  0.0e+00
    ## ma2      0.27070    0.000574    471.4391  0.0e+00
    ## omega   -1.79325    0.233210     -7.6894  0.0e+00
    ## alpha1   0.14348    0.030588      4.6906  3.0e-06
    ## beta1    0.35819    0.082956      4.3178  1.6e-05
    ## gamma1   0.41914    0.046728      8.9698  0.0e+00
    ## skew     1.32266    0.037586     35.1902  0.0e+00
    ## shape    3.54346    0.238225     14.8744  0.0e+00
    ## 
    ## LogLikelihood : 347.9765 
    ## 
    ## Information Criteria
    ## ------------------------------------
    ##                      
    ## Akaike       -0.22456
    ## Bayes        -0.20664
    ## Shibata      -0.22458
    ## Hannan-Quinn -0.21812
    ## 
    ## Weighted Ljung-Box Test on Standardized Residuals
    ## ------------------------------------
    ##                          statistic   p-value
    ## Lag[1]                      0.5812 4.459e-01
    ## Lag[2*(p+q)+(p+q)-1][8]     8.5925 3.969e-08
    ## Lag[4*(p+q)+(p+q)-1][14]   14.1511 4.171e-03
    ## d.o.f=3
    ## H0 : No serial correlation
    ## 
    ## Weighted Ljung-Box Test on Standardized Squared Residuals
    ## ------------------------------------
    ##                         statistic p-value
    ## Lag[1]                     0.4995  0.4797
    ## Lag[2*(p+q)+(p+q)-1][5]    1.1855  0.8164
    ## Lag[4*(p+q)+(p+q)-1][9]    2.4090  0.8510
    ## d.o.f=2
    ## 
    ## Weighted ARCH LM Tests
    ## ------------------------------------
    ##             Statistic Shape Scale P-Value
    ## ARCH Lag[3]    0.4215 0.500 2.000  0.5162
    ## ARCH Lag[5]    0.5974 1.440 1.667  0.8545
    ## ARCH Lag[7]    1.2835 2.315 1.543  0.8636
    ## 
    ## Nyblom stability test
    ## ------------------------------------
    ## Joint Statistic:  5.2333
    ## Individual Statistics:              
    ## ar1    0.63051
    ## ma1    1.18685
    ## ma2    1.11562
    ## omega  2.10211
    ## alpha1 0.08261
    ## beta1  2.07607
    ## gamma1 0.15883
    ## skew   0.33181
    ## shape  2.56140
    ## 
    ## Asymptotic Critical Values (10% 5% 1%)
    ## Joint Statistic:          2.1 2.32 2.82
    ## Individual Statistic:     0.35 0.47 0.75
    ## 
    ## Sign Bias Test
    ## ------------------------------------
    ##                    t-value    prob sig
    ## Sign Bias            1.600 0.10965    
    ## Negative Sign Bias   0.602 0.54725    
    ## Positive Sign Bias   2.540 0.01115  **
    ## Joint Effect         6.815 0.07804   *
    ## 
    ## 
    ## Adjusted Pearson Goodness-of-Fit Test:
    ## ------------------------------------
    ##   group statistic p-value(g-1)
    ## 1    20     20.37       0.3726
    ## 2    30     36.82       0.1510
    ## 3    40     45.07       0.2328
    ## 4    50     52.03       0.3567
    ## 
    ## 
    ## Elapsed time : 1.364722
    

    All coefficients are statistically significative. However, based on above reported Weighted Ljung-Box Test on Standardized Residuals p-values, we reject the null hypothesis of no correlation of residuals for the present model. Hence model ARMA(1,2)+eGARCH(1,1) is not able to capture all the structure of our time series.

    ARMA-GARCH: ARMA(2,3) + eGARCH(1,1)

    garchspec 
    ## 
    ## *---------------------------------*
    ## *          GARCH Model Fit        *
    ## *---------------------------------*
    ## 
    ## Conditional Variance Dynamics    
    ## -----------------------------------
    ## GARCH Model  : eGARCH(1,1)
    ## Mean Model   : ARFIMA(2,0,3)
    ## Distribution : sstd 
    ## 
    ## Optimal Parameters
    ## ------------------------------------
    ##         Estimate  Std. Error   t value Pr(>|t|)
    ## ar1     -0.18607    0.008580  -21.6873  0.0e+00
    ## ar2      0.59559    0.004596  129.5884  0.0e+00
    ## ma1     -0.35619    0.013512  -26.3608  0.0e+00
    ## ma2     -0.83010    0.004689 -177.0331  0.0e+00
    ## ma3      0.26277    0.007285   36.0678  0.0e+00
    ## omega   -1.92262    0.226738   -8.4795  0.0e+00
    ## alpha1   0.14382    0.033920    4.2401  2.2e-05
    ## beta1    0.31060    0.079441    3.9098  9.2e-05
    ## gamma1   0.43137    0.043016   10.0281  0.0e+00
    ## skew     1.32282    0.031382   42.1523  0.0e+00
    ## shape    3.48939    0.220787   15.8043  0.0e+00
    ## 
    ## Robust Standard Errors:
    ##         Estimate  Std. Error   t value Pr(>|t|)
    ## ar1     -0.18607    0.023940   -7.7724 0.000000
    ## ar2      0.59559    0.022231   26.7906 0.000000
    ## ma1     -0.35619    0.024244  -14.6918 0.000000
    ## ma2     -0.83010    0.004831 -171.8373 0.000000
    ## ma3      0.26277    0.030750    8.5453 0.000000
    ## omega   -1.92262    0.266462   -7.2154 0.000000
    ## alpha1   0.14382    0.032511    4.4239 0.000010
    ## beta1    0.31060    0.095329    3.2582 0.001121
    ## gamma1   0.43137    0.047092    9.1602 0.000000
    ## skew     1.32282    0.037663   35.1225 0.000000
    ## shape    3.48939    0.223470   15.6146 0.000000
    ## 
    ## LogLikelihood : 356.4994 
    ## 
    ## Information Criteria
    ## ------------------------------------
    ##                      
    ## Akaike       -0.22888
    ## Bayes        -0.20698
    ## Shibata      -0.22891
    ## Hannan-Quinn -0.22101
    ## 
    ## Weighted Ljung-Box Test on Standardized Residuals
    ## ------------------------------------
    ##                          statistic p-value
    ## Lag[1]                      0.7678 0.38091
    ## Lag[2*(p+q)+(p+q)-1][14]    7.7336 0.33963
    ## Lag[4*(p+q)+(p+q)-1][24]   17.1601 0.04972
    ## d.o.f=5
    ## H0 : No serial correlation
    ## 
    ## Weighted Ljung-Box Test on Standardized Squared Residuals
    ## ------------------------------------
    ##                         statistic p-value
    ## Lag[1]                      0.526  0.4683
    ## Lag[2*(p+q)+(p+q)-1][5]     1.677  0.6965
    ## Lag[4*(p+q)+(p+q)-1][9]     2.954  0.7666
    ## d.o.f=2
    ## 
    ## Weighted ARCH LM Tests
    ## ------------------------------------
    ##             Statistic Shape Scale P-Value
    ## ARCH Lag[3]     1.095 0.500 2.000  0.2955
    ## ARCH Lag[5]     1.281 1.440 1.667  0.6519
    ## ARCH Lag[7]     1.940 2.315 1.543  0.7301
    ## 
    ## Nyblom stability test
    ## ------------------------------------
    ## Joint Statistic:  5.3764
    ## Individual Statistics:              
    ## ar1    0.12923
    ## ar2    0.20878
    ## ma1    1.15005
    ## ma2    1.15356
    ## ma3    0.97487
    ## omega  2.04688
    ## alpha1 0.09695
    ## beta1  2.01026
    ## gamma1 0.18039
    ## skew   0.38131
    ## shape  2.40996
    ## 
    ## Asymptotic Critical Values (10% 5% 1%)
    ## Joint Statistic:          2.49 2.75 3.27
    ## Individual Statistic:     0.35 0.47 0.75
    ## 
    ## Sign Bias Test
    ## ------------------------------------
    ##                    t-value    prob sig
    ## Sign Bias           1.4929 0.13556    
    ## Negative Sign Bias  0.6317 0.52766    
    ## Positive Sign Bias  2.4505 0.01432  **
    ## Joint Effect        6.4063 0.09343   *
    ## 
    ## 
    ## Adjusted Pearson Goodness-of-Fit Test:
    ## ------------------------------------
    ##   group statistic p-value(g-1)
    ## 1    20     17.92       0.5278
    ## 2    30     33.99       0.2395
    ## 3    40     44.92       0.2378
    ## 4    50     50.28       0.4226
    ## 
    ## 
    ## Elapsed time : 1.660402
    

    All coefficients are statistically significative. No correlation within standardized residuals or standardized squared residuals is found. All ARCH effects are properly captured by the model. The Adjusted Pearson Goodness-of-fit test does not reject the null hypothesis that the empirical distribution of the standardized residuals and the selected theoretical distribution are the same. However:

    * the Nyblom stability test null hypothesis that the model parameters are constant over time is rejected for some of them

    * the Positive Sign Bias null hypothesis is rejected at 5% of significance level; this kind of test focuses on the effect of large and small positive shocks

    See ref. [11] for an explanation about GARCH model diagnostic.

    Some diagnostic plots are shown.

    par(mfrow=c(2,2))
    plot(garchfit, which=8)
    plot(garchfit, which=9)
    plot(garchfit, which=10)
    plot(garchfit, which=11)
    

    We show the original DJIA daily trade volume log ratio time series with the mean model fit (red line) and the conditional volatility (blue line).

    par(mfrow=c(1,1))
    cond_volatility 

    Model Equation

    Combining both ARMA(2,3) and eGARCH(1,1) model equations we have:

    [
    begin{equation}
    begin{cases}
    y_{t} – phi_{1}y_{t-1} – phi_{2}y_{t-2} = phi_{0} + u_{t} + theta_{1}u_{t-1} + theta_{2}u_{t-2} + theta_{3}u_{t-3}

    u_{t} = sigma_{t}epsilon_{t}, epsilon_{t}=N(0,1)

    ln(sigma_{t}^2) = omega + sum_{j=1}^{q} (alpha_{j} epsilon_{t-j}^2 + gamma (epsilon_{t-j} – E|epsilon_{t-j}|)) + sum_{i=1}^{p} beta_{i} ln(sigma_{t-1}^2)
    end{cases}
    end{equation}
    ]

    Using the model resulting coefficients we obtain:

    [
    begin{equation}
    begin{cases}
    y_{t} + 0.186 y_{t-1} – 0.595 y_{t-2} = u_{t} – 0.356u_{t-1} – 0.830 u_{t-2} + 0.263 u_{t-3}

    u_{t} = sigma_{t}epsilon_{t}, epsilon_{t}=N(0,1)

    ln(sigma_{t}^2) = -1.922 +0.144 epsilon_{t-1}^2 + 0.431 (epsilon_{t-1} – E|epsilon_{t-1}|)) + 0.310 ln(sigma_{t-1}^2)
    end{cases}
    end{equation}
    ]

    Daily Volume Log Ratio Volatility Analysis

    Here is the plot of conditional volatility as resulting from our model ARMA(2,2)+eGARCH(1,1).

    plot(cond_volatility)
    

    Line plots of conditional volatility by year are shown.

    par(mfrow=c(6,2))
    pl 

    Conditional volatility box-plots by year are shown.

    par(mfrow=c(1,1))
    cond_volatility_df 

    If you have any questions, please feel free to comment below.

    References

    1. Dow Jones Industrial Average
    2. Skewness ]
    3. Kurtosis ]
    4. An introduction to the analysis of financial data with R, Wiley, Ruey S. Tsay
    5. Time series analysis and its applications, Springer ed., R.H. Shumway, D.S. Stoffer
    6. Applied Econometric Time Series, Wiley, W. Enders, 4th ed.
    7. Forecasting – Principle and Practice, Texts, R.J. Hyndman
    8. Options, Futures and other Derivatives, Pearson ed., J.C. Hull
    9. An introduction to rugarch package
    10. Applied Econometrics with R, Achim Zeileis, Christian Kleiber – Springer Ed.
    11. GARCH modeling: diagnostic tests

    Disclaimer

    Any securities or databases referred in this post are solely for illustration purposes, and under no regard should the findings presented here be interpreted as investment advice or promotion of any particular security or source.

    Related Post

    1. Dow Jones Stock Market Index (3/4): Log Returns GARCH Model
    2. Leaf Plant Classification: Statistical Learning Model – Part 2
    3. NYC buses: C5.0 classification with R; more than 20 minute delay?
    4. NYC buses: Cubist regression with more predictors
    5. NYC buses: simple Cubist regression

    To leave a comment for the author, please follow the link and comment on their blog: R Programming – DataScience+.

    R-bloggers.com offers daily e-mail updates about R news and tutorials on topics such as: Data science, Big Data, R jobs, visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more…

    Comments are closed.