# additional nonlinear terms in the EGARCH model, maximum likelihood estimation

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.

• It would help to comment on this if you made your data and code available. However, in general, my advice would be to fit the model in a number of phases. IN ths first phase you should just fit the model with the parameters which are responsible for the main variation in the variance. Also the $w_i$ part of the equation doesn't look right. Finally use $\ln$ and not $ln$. Math functions should not use math italics, it offends the eye. Commented Dec 11, 2016 at 15:53
• Looks much better. Still no data or equations. I can knock together one of these models in 15 minutes with software that uses automatic differentiation to get exact derivatives for all the model parameters if I know exactly what the equations are, and run it if I have the data. Lacking that one is just wasting time. Commented Dec 13, 2016 at 18:40
• Thanks Dave!! The dataset is rather big. Can I sent it you via e-mail? I have used Eviews for coding. I noticed you use R. So I'm not sure if its helping to post my tedious code here.. Commented Dec 14, 2016 at 15:52
• email is good. my junk email (so that I can post the address) is [email protected] Commented Dec 14, 2016 at 16:09

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 = 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.

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

• Thanks Dave!! I still have a couple of questions. 1. You worte "the model almost certainly has too many parameters". But we can still somehow fit the model. So in this case what should we do? Do you suggest to keep the model or drop some of the the $S_{j,t}$ variables or find another process for the variance? Commented Dec 16, 2016 at 14:27
• 2. Is it okay just to restrict the $\omega_{3}$ parameter to be strictly smaller than one? I mean, if we know that the conditional variance might be fractionally integrated, wouldn't we make type II error by restricting $\omega_{3}<1$? 3. If it is safe to restrict $\omega_{3}$, can I use the estimated process to perform further calculations? For example counterfacual analysis? Commented Dec 17, 2016 at 12:25
• Well it turns out that starting from the solution where $\omega_3$ is constrained to be < 0.98 and relaxing this constraint to be < 1.02 the model converges to a solution with $\omega_3=0.995013839900$. So I think that solves the problem. I'll put the solution in the answer. Commented Dec 17, 2016 at 16:11