2

I am trying to do the restriction test for GARCH model (ugarch from 'rugarch' package) using the following hypothesis:

 H0: alpha1 + beta1 = 1

 H1: alpha1 + beta1 ≠ 1 

So I am trying to follow the advice from https://stats.stackexchange.com/questions/151573/testing-the-sum-of-garch1-1-parameters/151578?noredirect=1#comment629951_151578

1.Specify the restricted model using ugarchspec with option variance.model = list(model = "sGARCH") and estimate it using ugarchfit. Obtain the log-likelihood from the slot fit sub-slot likelihood.

2.Specify the restricted model using ugarchspec with option variance.model = list(model = "iGARCH") and estimate it using ugarchfit. Obtain the log-likelihood as above.

3.Calculate LR=2(Log-likelihood of unrestricted model − Log-likelihood of restricted model) and Obtain the p-value as pchisq(q = LR, df = 1).

I have the following 'sGARCH' and 'iGARCH' models I am using from 'rugarch' package.

(A) sGARCH (unrestricted model):

 speccR = ugarchspec(variance.model=list(model = "sGARCH",garchOrder=c(1,1)),mean.model=list(armaOrder=c(0,0), include.mean=TRUE,archm = TRUE, archpow = 1))

 ugarchfit(speccR, data=data.matrix(P),fit.control = list(scale = 1))

And the following is this sGARCH output:

    *---------------------------------*
    *          GARCH Model Fit        *
    *---------------------------------*

    Conditional Variance Dynamics   
    -----------------------------------
    GARCH Model     : sGARCH(1,1)
    Mean Model      : ARFIMA(0,0,0)
    Distribution    : norm 

    Optimal Parameters
    ------------------------------------
            Estimate  Std. Error  t value Pr(>|t|)
    mu     -0.000355    0.001004 -0.35377 0.723508
    archm   0.096364    0.039646  2.43059 0.015074
    omega   0.000049    0.000010  4.91096 0.000001
    alpha1  0.289964    0.021866 13.26117 0.000000
    beta1   0.709036    0.023200 30.56156 0.000000

    Robust Standard Errors:
            Estimate  Std. Error  t value Pr(>|t|)
    mu     -0.000355    0.001580 -0.22482 0.822122
    archm   0.096364    0.056352  1.71002 0.087262
    omega   0.000049    0.000051  0.96346 0.335316
    alpha1  0.289964    0.078078  3.71375 0.000204
    beta1   0.709036    0.111629  6.35173 0.000000

    LogLikelihood : 5411.828 

    Information Criteria
    ------------------------------------

    Akaike       -3.9180
    Bayes        -3.9073
    Shibata      -3.9180
    Hannan-Quinn -3.9141

    Weighted Ljung-Box Test on Standardized Residuals
    ------------------------------------
                            statistic p-value
    Lag[1]                      233.2       0
    Lag[2*(p+q)+(p+q)-1][2]     239.1       0
    Lag[4*(p+q)+(p+q)-1][5]     247.4       0
    d.o.f=0
    H0 : No serial correlation

    Weighted Ljung-Box Test on Standardized Squared Residuals
    ------------------------------------
                            statistic p-value
    Lag[1]                      4.695 0.03025
    Lag[2*(p+q)+(p+q)-1][5]     5.941 0.09286
    Lag[4*(p+q)+(p+q)-1][9]     7.865 0.13694
    d.o.f=2

    Weighted ARCH LM Tests
    ------------------------------------
                Statistic Shape Scale P-Value
    ARCH Lag[3]     0.556 0.500 2.000  0.4559
    ARCH Lag[5]     1.911 1.440 1.667  0.4914
    ARCH Lag[7]     3.532 2.315 1.543  0.4190

    Nyblom stability test
    ------------------------------------
    Joint Statistic:  5.5144
    Individual Statistics:             
    mu     0.5318
    archm  0.4451
    omega  1.3455
    alpha1 4.1443
    beta1  2.2202

    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.2384 0.8116    
    Negative Sign Bias  1.1799 0.2381    
    Positive Sign Bias  1.1992 0.2305    
    Joint Effect        2.9540 0.3988    


    Adjusted Pearson Goodness-of-Fit Test:
    ------------------------------------
      group statistic p-value(g-1)
    1    20     272.1    9.968e-47
    2    30     296.9    3.281e-46
    3    40     313.3    1.529e-44
    4    50     337.4    1.091e-44


    Elapsed time : 0.4910491 

(B) iGARCH (restricted model):

 speccRR = ugarchspec(variance.model=list(model = "iGARCH",garchOrder=c(1,1)),mean.model=list(armaOrder=c(0,0), include.mean=TRUE,archm = TRUE, archpow = 1))

 ugarchfit(speccRR, data=data.matrix(P),fit.control = list(scale = 1))

However, I get the following output of beta1 with N/A in its standard error, t-value, and p-value.

The following is the iGARCH output:

    *---------------------------------*
    *          GARCH Model Fit        *
    *---------------------------------*

    Conditional Variance Dynamics   
    -----------------------------------
    GARCH Model     : iGARCH(1,1)
    Mean Model      : ARFIMA(0,0,0)
    Distribution    : norm 

    Optimal Parameters
    ------------------------------------
            Estimate  Std. Error  t value Pr(>|t|)
    mu     -0.000355    0.001001 -0.35485 0.722700
    archm   0.096303    0.039514  2.43718 0.014802
    omega   0.000049    0.000008  6.42826 0.000000
    alpha1  0.290304    0.021314 13.62022 0.000000
    beta1   0.709696          NA       NA       NA

    Robust Standard Errors:
            Estimate  Std. Error  t value Pr(>|t|)
    mu     -0.000355    0.001488  -0.2386 0.811415
    archm   0.096303    0.054471   1.7680 0.077066
    omega   0.000049    0.000032   1.5133 0.130215
    alpha1  0.290304    0.091133   3.1855 0.001445
    beta1   0.709696          NA       NA       NA

    LogLikelihood : 5412.268 

    Information Criteria
    ------------------------------------

    Akaike       -3.9190
    Bayes        -3.9105
    Shibata      -3.9190
    Hannan-Quinn -3.9159

    Weighted Ljung-Box Test on Standardized Residuals
    ------------------------------------
                            statistic p-value
    Lag[1]                      233.2       0
    Lag[2*(p+q)+(p+q)-1][2]     239.1       0
    Lag[4*(p+q)+(p+q)-1][5]     247.5       0
    d.o.f=0
    H0 : No serial correlation

    Weighted Ljung-Box Test on Standardized Squared Residuals
    ------------------------------------
                            statistic p-value
    Lag[1]                      4.674 0.03063
    Lag[2*(p+q)+(p+q)-1][5]     5.926 0.09364
    Lag[4*(p+q)+(p+q)-1][9]     7.860 0.13725
    d.o.f=2

    Weighted ARCH LM Tests
    ------------------------------------
                Statistic Shape Scale P-Value
    ARCH Lag[3]    0.5613 0.500 2.000  0.4538
    ARCH Lag[5]    1.9248 1.440 1.667  0.4881
    ARCH Lag[7]    3.5532 2.315 1.543  0.4156

    Nyblom stability test
    ------------------------------------
    Joint Statistic:  1.8138
    Individual Statistics:             
    mu     0.5301
    archm  0.4444
    omega  1.3355
    alpha1 0.4610

    Asymptotic Critical Values (10% 5% 1%)
    Joint Statistic:         1.07 1.24 1.6
    Individual Statistic:    0.35 0.47 0.75

    Sign Bias Test
    ------------------------------------
                       t-value   prob sig
    Sign Bias           0.2252 0.8218    
    Negative Sign Bias  1.1672 0.2432    
    Positive Sign Bias  1.1966 0.2316    
    Joint Effect        2.9091 0.4059    


    Adjusted Pearson Goodness-of-Fit Test:
    ------------------------------------
      group statistic p-value(g-1)
    1    20     273.4    5.443e-47
    2    30     300.4    6.873e-47
    3    40     313.7    1.312e-44
    4    50     337.0    1.275e-44


    Elapsed time : 0.365 

If I calculate the log-likelihood difference to derive the chi-square value as suggested I get negative value as such:

 2*(5411.828-5412.268)=-0.88

The Log-likelihood of the restricted model "iGARCH" is 5412.268 which is higher than the Log-likelihood of the unrestricted model "sGARCH" of 5411.828 which should not happen.

The data I use are as follows in time series manner (I am only posting first 100 data due to space limit):

   Time      P
    1   0.454213593
    2   0.10713195
    3   -0.106819399
    4   -0.101610699
    5   -0.094327846
    6   -0.037176107
    7   -0.101550977
    8   -0.016309894
    9   -0.041889484
    10  0.103384357
    11  -0.011746377
    12  0.063304432
    13  0.059539249
    14  -0.049946177
    15  -0.023251656
    16  0.013989353
    17  -0.002815588
    18  -0.009678745
    19  -0.011139779
    20  0.031592303
    21  -0.02348106
    22  -0.007206591
    23  0.077422089
    24  0.064632768
    25  -0.003396734
    26  -0.025524166
    27  -0.026632474
    28  0.014614485
    29  -0.012380888
    30  -0.007463018
    31  0.022759969
    32  0.038667465
    33  -0.028619484
    34  -0.021995984
    35  -0.006162809
    36  -0.031187399
    37  0.022455611
    38  0.011419264
    39  -0.005700445
    40  -0.010106343
    41  -0.004310162
    42  0.00513715
    43  -0.00498106
    44  -0.021382251
    45  -0.000694252
    46  -0.033326085
    47  0.002596086
    48  0.011008057
    49  -0.004754233
    50  0.008969559
    51  -0.00354088
    52  -0.007213115
    53  -0.003101495
    54  0.005016228
    55  -0.010762641
    56  0.030770993
    57  -0.015636325
    58  0.000875417
    59  0.03975863
    60  -0.050207219
    61  0.011308261
    62  -0.021453315
    63  -0.003309127
    64  0.025687191
    65  0.009467306
    66  0.005519485
    67  -0.011473758
    68  0.00223934
    69  -0.000913651
    70  -0.003055385
    71  0.000974694
    72  0.000288611
    73  -0.002432251
    74  -0.0016975
    75  -0.001565034
    76  0.003332848
    77  -0.008007295
    78  -0.003086435
    79  -0.00160435
    80  0.005825885
    81  0.020078093
    82  0.018055453
    83  0.181098137
    84  0.102698818
    85  0.128786594
    86  -0.013587077
    87  -0.038429879
    88  0.043637258
    89  0.042741709
    90  0.016384872
    91  0.000216317
    92  0.009275681
    93  -0.008595197
    94  -0.016323335
    95  -0.024083247
    96  0.035922206
    97  0.034863621
    98  0.032401779
    99  0.126333922
    100 0.054751935

In order to perform the restriction test from my H0 and H1 hypothesis, may I know how I can fix this problem?

Eric
  • 528
  • 1
  • 8
  • 26
  • Does anyone know the answer especially when I am doing restriction test using GARCH-M models? – Eric Mar 12 '18 at 01:04

2 Answers2

1

There seems to be a problem with the estimation procedure... Since one model is a restricted version of the other, using iGARCH should indeed lead to a lower likelihood.

Using the subset of your data,

fit1 <- ugarchfit(speca, data = data.matrix(P)) 
# [1] 161.7373
fit2 <- ugarchfit(speca2, data = data.matrix(P))
# [1] 165.333

As I said in my deleted post, those numbers looked suspicious, as if they are -loglikelihoods. However, recovering the likelihood from the residuals gives

-sum(log(2 * pi * sigma(fit1)^2)) / 2 - sum(residuals(fit1, standardize = TRUE)^2) / 2
# [1] 161.7373
-sum(log(2 * pi * sigma(fit2)^2)) / 2 - sum(residuals(fit2, standardize = TRUE)^2) / 2
# [1] 165.333

Meaning that my suspicion was wrong (it must then be that the density values are > 1). For this reason, I think there is no way to use the current output to construct a test. The iGARCH restriction fits miraculously well..

However, some experimenting showed that using

fit.control = list(scale = 1)

changes things. In particular,

fit1 <- ugarchfit(speca, data = data.matrix(P), fit.control = list(scale = 1))
likelihood(fit1)
# [1] 161.7373
-sum(log(2 * pi * sigma(fit1)^2)) / 2 - sum(residuals(fit1, standardize = TRUE)^2) / 2
# [1] 161.7373

fit2 <- ugarchfit(speca2, data = data.matrix(P), fit.control = list(scale = 1))
likelihood(fit2)
# [1] 19.5233
-sum(log(2 * pi * sigma(fit2)^2)) / 2 - sum(residuals(fit2, standardize = TRUE)^2) / 2
# [1] 19.5233

That would somewhat make sense given

(page 25) "scaling sometimes facilitates the estimation process"

(page 46) Q: My model does not converge, what can I do?

"...Additionally, in the fit.control list of the fitting routines, the option to perform scaling of the data prior to fitting usually helps, although it is not available under some setups..."

However, it again is suspicious that the likelihood of the first model remains the same. Then we have that

fit1 <- ugarchfit(speca, data = data.matrix(P), fit.control = list(scale = 0), solver.control = list(trace = TRUE))
# 
# Iter: 1 fn: -161.7373  Pars:  -0.0454619  0.0085993  0.0002706  0.0593231  # 0.6898473
# Iter: 2 fn: -161.7373  Pars:  -0.0454619  0.0085993  0.0002706  0.0593231  0.6898473
# solnp--> Completed in 2 iterations
coef(fit1)
#            mu        mxreg1         omega        alpha1         beta1 
# -0.0454619274  0.0085992743  0.0002706018  0.0593231138  0.6898472858 
fit1 <- ugarchfit(speca, data = data.matrix(P), fit.control = list(scale = 1), solver.control = list(trace = TRUE))

# Iter: 1 fn: 114.8143   Pars:  -0.72230  0.13663  0.06830  0.05930  0.68988
# Iter: 2 fn: 114.8143   Pars:  -0.72228  0.13662  0.06830  0.05931  0.68986
# solnp--> Completed in 2 iterations
coef(fit1)
#           mu       mxreg1        omega       alpha1        beta1 
# -0.045463099  0.008599494  0.000270610  0.059310622  0.689858216

and

fit2 <- ugarchfit(speca2, data = data.matrix(P), fit.control = list(scale = 0), solver.control = list(trace = TRUE))

# Iter: 1 fn: -165.3330  Pars:   0.0292439 -0.0051098  0.0002221  0.7495846
# Iter: 2 fn: -165.3330  Pars:   0.0292434 -0.0051097  0.0002221  0.7495853
# solnp--> Completed in 2 iterations
coef(fit2)
#            mu        mxreg1         omega        alpha1         beta1 
#  0.0292434276 -0.0051096984  0.0002221457  0.7495853224  0.2504146776 
fit2 <- ugarchfit(speca2, data = data.matrix(P), fit.control = list(scale = 1), solver.control = list(trace = TRUE))

# Iter: 1 fn: 111.2185   Pars:   0.46462 -0.08118  0.05607  0.74959
# Iter: 2 fn: 111.2185   Pars:   0.46458 -0.08118  0.05607  0.74959
# solnp--> Completed in 2 iterations
coef(fit2)
#          mu      mxreg1       omega      alpha1       beta1 
#  0.46458110 -0.08117626  0.05607215  0.74959242  0.25040758 

Which makes things even stranger due to multiple inconsistencies...

Julius Vainora
  • 47,421
  • 9
  • 90
  • 102
  • Is there a way that I can post the full data set I am using for you to replicate my code and result? I want to attach a notepad file but cannot. – Eric Mar 16 '18 at 22:15
  • @Eric, I don't think you can do that directly; two alternatives that come to mind are pastebin.com and dropbox. – Julius Vainora Mar 16 '18 at 23:17
0

This is the answer I received from the package author "Alexios Galanos":

The problem is that there is a restriction on the stationarity of the GARCH model which may interfere with the solver convergence for models which are on the border of stationarity. Here is the solution:

  library(rugarch)
  library(xts)
  dat<-read.table("data.txt",header = TRUE, stringsAsFactors = FALSE)
  dat = xts(dat[,2], as.Date(strptime(dat[,1],"%d/%m/%Y")))

  spec1<-ugarchspec(mean.model=list(armaOrder=c(0,0), archm=TRUE, archpow=1), variance.model=list(model="iGARCH"))
  spec2<-ugarchspec(mean.model=list(armaOrder=c(0,0), archm=TRUE, archpow=1), variance.model=list(model="sGARCH"))
  mod1<-ugarchfit(spec1, dat, solver="solnp")
  mod2<-ugarchfit(spec2,dat)
  persistence(mod2)
  >0.999

 # at the limit of the internal constraint

 mod2<-ugarchfit(spec2, dat, solver="solnp", fit.control = list(stationarity=0))
  likelihood(mod2)
  >5428.871


  likelihood(mod1)

  >5412.268
  persistence(mod2)
  1.08693
  # above the limit

  Here is one solution to change the constraint:

  .garchconbounds2= function(){
    return(list(LB = 1e-12,UB = 0.99999999999))
  }
  assignInNamespace(x = ".garchconbounds", value=.garchconbounds2, ns="rugarch")
  mod2<-ugarchfit(spec2, dat, solver="solnp")

  likelihood(mod2)
  >5412.268

Now the value is the same as the constrained model (they are both effectively integrated), but the constrained model has one less parameter to estimate.

I don't even need a fit.control=list(scale=1) at all here. Probably better to delete this scale.

Eric
  • 528
  • 1
  • 8
  • 26