additional nonlinear terms in the EGARCH model, maximum likelihood estimation

Question

I have estimated the an EGARCH model of the following form via maximum likelihood.

$\Delta y_{t} = \alpha + \sum^{12}_{j}\beta_{j}S_{j,t}(1 + \gamma_{1}D_{1,t} + \gamma_{2}D_{2,t}) + \varepsilon_{t}$

$\ln(\sigma^2_{t}) = \omega_{0} + \sum^{2}_{i=1} \omega_{i} |\frac{\varepsilon_{t-i}}{\sigma_{t-i}}| + \lambda \frac{\varepsilon_{t-1}}{\sigma_{t-1}} + \omega_{3} \ln(\sigma^{2}_{t-1}) + \sum^{12}_{j}\delta_{j}|S_{j,t}|(1 + \rho_{1}D_{1,t} + \rho_{2}D_{2,t})$

where $D_{1}$ and $D_{2}$ are dummies, $\Delta y_{t}$ is the change in dollar-euro spot rates and $t$ indexes days.

Everytime I get different results after estimating the model. My guess is maximum likelihood estimation performs poorly at handling so many nonlinear terms. But I cannot find a theoretical justification for this. Can anyone give me some hint? Looking forward for discussions.

Answer 1

A bit of philosophy about model fitting. I used to get into arguments with the R cognoscenti, some of whom maintained that the failure of an unstable algorithm to fit the model should be celebrated as a useful diagnostic for detecting uninformative data or poor model specification. They never explained why this usefulness of failure should not be extended to linear models for which there exist extremely stable algorithms that can fit almost any model. The idea is that with linear models you first fit the model and then rely on the fitted model produce diagnostics which can inform you about the questions of reliability. I think the same approach should be employed for nonlinear models. It is just a lot more difficult.

This model almost certainly has too many parameters, but lets see if it can be fitted with all of them in some reasonable way. The first problem is to deal with the initial values for the log-variance. One way which I chose is to pick the equilibrium condition assuming the $\sigma^2_{t-1}=\sigma^2_t$ so that

$$\ln(\sigma_t^2) = \omega_{0} + \omega_{3} \ln(\sigma^{2}_{t})$$

so that $$\ln(\sigma^2_1)=\frac{\omega_0}{1-\omega_3}$$

and then you see the first problem is that for the model to be stationary one needs the condition $\omega_3<1$. I have bounded it to be < 0.98.

This assumption seems to worry the OP, so I will also investigate relaxing this condition.

With that condition it appears to converge quite nicely.

Now one can fit the model leaving out various terms and doing a likelihood ratio test. for example these are the results with and without $\omega_2$

Number of parameters = 33 Objective function value = -63.4744 Maximum gradient component = 0.000247846

Number of parameters = 34 Objective function value = -64.4108 Maximum gradient component = 7.15130e-05

So the log-likelihood difference is 64.4108 - 63.4744 which is not significant.

Ths is the AD MOdel Builder code for the model.

DATA_SECTION
  init_int n
  init_int m
  init_matrix data(1,n,1,m)
  matrix S(1,n,1,12)
  vector deltaE
  matrix D(1,n,1,2)
LOC_CALCS
  for (int i=1;i<=n;i++)
  {
    S(i)=data(i)(5,16).shift(1);
    D(i)=data(i)(17,18).shift(1);
  }
  deltaE=column(data,3);
PARAMETER_SECTION
  init_number alpha
  init_vector beta(1,12,2)
  init_vector gamma(1,2,3)
  init_number omega0
  init_vector omega(1,2,3)
  init_number lambda(3)
  init_bounded_number omega3(0.01,0.98,2)
  init_vector delta(1,12,4)
  init_vector rho(1,2,5)
  vector s(1,n)
  vector logsigma2(1,n)
  vector eta(1,n)
  vector pred_deltaE(1,n)
  objective_function_value f
PROCEDURE_SECTION
  double ee=1.e-5;
  if (current_phase()>3)
   ee=1.e-15;
  logsigma2(1)=omega0/(1.0-omega3);
  logsigma2(2)=omega0/(1.0-omega3);
  s(1)=mfexp(0.5*logsigma2(1));
  s(2)=mfexp(0.5*logsigma2(2));
  pred_deltaE(1)=alpha+beta*S(1)*(1.0+gamma*D(1));
  eta(1)=deltaE(1)-pred_deltaE(1);
  for (int i=2;i<=n;i++)
  {
    pred_deltaE(i)=alpha+beta*S(i)*(1.0+gamma*D(i));
    if (i>2)
      logsigma2(i)=omega0
       +omega(1)*sfabs(eta(i-1)/(ee+s(i-1)))
       +lambda*eta(i-1)/(ee+s(i-1))
       +omega3*logsigma2(i-1)
       +delta*fabs(S(i))*(1+rho*D(i))
       +omega(2)*sfabs(eta(i-2)/(ee+s(i-2)));

   s(i)=mfexp(0.5*logsigma2(i));  
   eta(i)=deltaE(i)-pred_deltaE(i);
 }

 f+=0.5*sum(logsigma2) +  0.5*sum(square(elem_div(eta,(ee+s))));  

REPORT_SECTION
  report << "s" << endl << s << endl;

The plot of the variance looks like

So it appears that the variances are too large at the beginning. One can try letting them be free parameters for periods 1 and 2. Unfortunately this produces even larger variances for the first periods.

So it is not difficult to fit this model, but one needs some real world experience with this kind of problem to interpret the results.

These are the parameter estimates for the model assuming the initial equilibrium condition.

# alpha:
 -0.0251881678332 
# beta:
 0.0209472699255 -0.0587572424810 0.0181180949961 -0.488987381552
-0.0234483021308 0.0263686779804 -0.0189694436772 -0.0228941063294
-0.168297324565 -0.0133886102171 -0.0859940692792 0.0213416534723
# gamma:
  -0.501910563068 1.46851632417
# omega0:
 0.00127933766942
# omega(1,2):
 -0.0914689003618 0.0741717590471
# lambda:
 -0.0248365784280
# omega3:
 0.979999999502
# delta:
 -0.0287602387730 0.00157907201125 0.0130661282781 -0.0860609108614
 0.00229332323984 0.00679611779048 -0.0153838157778 0.0203405212533
 -0.00284365063170 0.00331155914496 0.0555029674112 -0.177277000904
# rho:
  2.15203087781 6.01624005387

Now lets investigate relaxing the condition that $\omega_3<1.0$

The easiest thing to do is to fix $\omega_3=1$ and see what happens. Since we can no longer use the equilibrium condition for $\sigma_1^2$, and $\sigma_2^2$ we will simply set them equal to a new parameter. So this add one more parameter to the model. What we find is that the model "blows up". So if one want to investigate non stationary solutions, a different approach to this maximum likelihood formulation will need to be employed. I think that is beyond the scope of the original question.

However one can get a sense of the flavour of the non stationary solutions by setting $\omega_3$ close to 1.0, say 0.995, if this converges to a solution. This model converges with no problem. A plot of the standard deviation is

Final Remarks

I don't think there is an easy "the solution" to this problem. My point is that it is possible to fit these models reliably using flexible software which makes it easy to change the model structure in different ways. Presumably in the real world one would use the time series of standard deviations to do something or other and could investigate whether the differences in the estimates from the various models are relevant.

Addendum:

It turns out that starting from the solution with $\omega_3<0.98$ one can find an unconstrained solution with $\omega_3=0.995013839900$ So I was a bit to pessimistic. The solution is

 # Number of parameters = 35  Objective function value = -89.3242  Maximum gradient component = 7.45641e-05
 # alpha:
 -0.0184773755820
 # beta:
  0.0776272250120 -0.0282086700728 0.0105870775566 -0.796727950111 -0.0329451864914 0.0574720791322 -0.00769006016788 0.0634996913172 -0.271844797702 -0.00148927676168 -0.107749917570 -0.0235315136652
 # gamma:
  -1.28418932915 0.253708104907
 # omega0:
 -0.0288892906469
 # omega:
  -0.0390792976332 0.0693366331859
 # lambda:
 -0.0174406066104
 # ls1:    log standard deviations for periods 1 and 2
 0.521640582874
 # omega3:
 0.995013839900
 # delta:
  -0.0454107838910 0.0442755011514 0.0848044432830 -0.202324778947 -0.0184050276271 0.00891208520744 0.00954099100492 0.00610579281349 -0.0181269442505 0.0502826182040 -0.0135121954742 -0.143376988696
 # rho:
  0.144999774670 1.88445900772

The plot of the standard deviations is