Financial Econometrics (EC7079), 2016
Computer class 4
Modelling time-variant conditional volatility
&
Simulating GARCH process in gretl
November 30, 2016
1 Introduction
Open the session file StockPrice.gretl. This contains data as well as a few results as a starting point
for further analysis.
The data set contains daily stock price on three companies: Apple, Microsoft, and HSBC plc.
They are named P ta, P tm, and P th respectively. All are trading on New York Stock Exchange,
and the prices are recorded in U.S. dollars.
Their log returns have been computed and included in your dataset. r ta, r tm, and r th denote
the log returns on Apple stocks, Microsoft stocks, and HSBC stocks respectively. In this class, you
will be exploring the properties of the log returns on Microsoft stocks. I shall leave it to you to
estimate a variety of models on the other two in your own time.
A plot of r tm, saved in your session file and printed below, shows that the series is certainly
stationary—it never seems to stray too far from 0. However, it also shows a typical property of
financial asset series: volatility clustering. There are calmer periods, in particular around 2005, but
there are also volatile periods; around year 2000, and the time of the immediate aftermath of the
Financial Crisis of 2008. There is a peculiar outlier early on in the sample range. The return on
Microsoft stock was down by 35% in one day on 19th of October 1987. This is the Black Monday
of October 1987, when stock markets around the world crashed. Two days later, the return on
Microsoft shares climbed back up by 16.5%.
The correlogram of r tm (not printed here but saved in your session file by the name of corr rt m)
indicate that the series might be generated by a mixed ARMA process, as neither ˆ ⇢k nor “ˆkk show
a clear-cut truncation. As an initial guess, ARMA(2,1) was estimated on the series. The estimation
result is printed on the next page (saved in your session file as an icon ARMA 21).
All parameter estimates are statistically significant. The residual correlogram is printed below
(saved in your session file as corr uhat1, where series uhat1 is the residual series obtained from
the estimation of ARMA(2,1) on r tm). Neither the residual correlogram nor its partial ACF
indicate any serious mis-specification. They indicate that the residuals are pretty much a white
noise process. However, the correlogram of the squares of residuals tells di↵erent story. I shall leave
it to you to create the squares of uhat1 and plot its correlogram. It shows that the variance of the
residuals is correlated across time.
1
2 Estimating ARCH & GARCH models
Although the residuals from the estimation of ARMA(2,1) may be a white noise, it appears from
the graph of r tm and the correlogram of the squared residuals that their volatility could well be
correlated. Testing for ARCH e↵ect should verify if the error term is generated by an ARCH
process.
2
2.1 Testing for ARCH errors
Recall from the lectures that the ARCH test developed by Engle (1982) is, in its simplest form, an
F-test on joint restrictions on the following regression model:
✏ˆ2
t = ↵0 + ↵1✏ˆ2 t!1 + ↵2✏ˆ2 t!2 + + · · · + ↵m✏ˆ2

(2)
The test statistic follows a %2 distribution asymptotically with m degrees of freedom if the null
hypothesis is true.
To conduct this test in gretl, double-click on the ARMA 21 estimation output icon in your icon
view window to open it. Click on Tests, then ARCH.
You should get ARCH test dialogue box, asking for the lag length in test equation (1). Leave it as
the default of 5, and click on OK to confirm. You should obtain the following result.
All 5 lags are highly significant. Even the fifth lag ˆ ✏2 t!5 has a p-value that is practically 0 and a
t-ratio of over 17. Needless to say, the %2 test emphatically rejects the null hypothesis.
3
2.2 Estimating ARCH model
The above test result indicates that we need to incorporate ARCH errors in any models we fit on
the log returns on Microsoft stocks. For now, we shall retain the two AR components as they were
highly significant from estimating ARMA(2,1). Instead of specifying MA(1) for the error term, we
shall let ✏t be generated by an ARCH process, without the MA component. This is partly because
there is no automatic clickable option in gretl for estimating ARMA(p, q) model with ARCH errors.
Let yt denote the log return r tm at time t. The model we wish to fit is defined as follows:
yt = c + “1yt!1 + “2yt!2 + ✏t (3)
where
✏t = phtvt; vt ⇠ iin(0, 1) (4)
ht = ↵0 +
mX i
=1
↵i✏2 t!i (5)
where the following parameter restrictions are applied:
↵0 > 0; ↵i ” 0, i = 1, 2, · · · , m; and
mX i
=1
↵i < 1
So the model is AR(2) with ARCH(m) errors. To estimate this model, on the main gretl window,
click on Model #! Time series #! GARCH…. Recall that ARCH(m) process is a restricted version
of GARCH(m, r)—the former is GARCH(m, 0). You will get a specify model dialogue box for
GARCH shown below.
4
Select r tm as the dependent variable. Click on lags at the bottom left of the dialogue box.
In the lag order dialogue box that follows, check the box for Lags of dependent variable, and select
order from 1 to 2.
Click on OK to return to the GARCH specify model dialogue box. You should see two lags
of r tm added to the box of Regressors. As for the GARCH orders, select 0 for GARCH p: and 4
for ARCH q:. Although the Lagrange Multiplier test above indicated that the fifth lag was highly
significant, 4 is the maximum lag gretl permits. This is because ARCH(m) with a log lag is very
well approximated by a GARCH process of small orders as we will see later on. Click on OK to
estimate the model. You should obtain the result printed below. Make sure to save this estimation
results by clicking on File, then save to session as an icon.
Notice that the statistical significance of the two AR parameters has all but disappeared once
we include ARCH errors (I shall leave it to you to think about why this is the case). As we would
expect, all four ARCH terms are highly significant. All three criteria select this AR(2)-ARCH(4)
model over ARMA(2,1) overwhelmingly—even Schwart-Bayesian Criterion, which always heavily
penalises ove-parameterisation, gives the ARCH model a much smaller number. Recall that these
criteria are loss functions (they are all varying functions of the magnitude of the residual variance),
therefore, the lower the number (greater in magnitude), the better the fit of the model.
Save the residuals from this model by clicking on Save, then Residuals.
In the variable attributes dialogue box, leave everything as the default with the name of uhat2, click
on OK. You should see a new series uhat2 added to your dataset. This model is the 2nd model as
far as gretl is concerned, because ARMA(2,1) has been estimated and saved in the session file.
5
gretl computes the estimate of the conditional variance series, ht, calling it ‘predicted error
variance’. In the window showing the AR(2)-ARCH(4) estimation result, click on Save, then
Predicted error variance:
gretl should open up variable attributes dialogue box, with the default name of h2 with description
fitted variance from model 2. Leave the name and description as the default and click on OK. You
should now see another series h2 added to your dataset. Double-click on h2, and in the data window
that follows, click on the graph icon to plot the series. Here is what you should obtain:
The changing volatility of the series noted in the introduction seems to have been well captured
by this model. Around 2005, the conditional volatility is very small compared to other time
periods. On the other hand, shortly after the year 2000, and after the Financial Crisis of 2008, we
see substantial increase in estimated ht, as we would expect.
2.3 Diagnostic checking
Although ARCH(4) seems inadequate from the LM test results, if most of the changing volatility
in ✏t has been captured by h˜t (estimated series of ht, h2), the error term which drives the entire
process, vt, should be a white noise. Recall from (4) that ✏t = phtvt. Then vt = ✏t/pht. We shall
recover vt using the estimates of ✏t and ht as below:
v˜t =
✏˜t
ph˜t
6
where ˜ ✏t denotes the residual series we have saved above as uhat2.
To create ˜ vt, click on Add, then Define new variable… to obtain the Add var dialogue box. Type
v2=uhat2/sqrt(h2):
and click on OK to confirm. Select the newly created series v2 in your main gretl window, click on
Variable, then Correlogram to obtain the following with the default number of lags of 38:
There is no indication of serial correlation in ˜ vt.
In your own time, create the square of residuals from ARMA(2,1) estimation on r tm (this
residual series is saved in your data file as uhat1) and examine the correlogram of ˆ ✏2 t . Also create
v˜2
t from the estimation of AR(2)-ARCH(4) and plot its correlogram. Compare the two. What do
you find? And what does the di↵erence between between the correlograms of ˆ ✏2 t and ˜ vt2 tell you?
2.4 Estimating GARCH model
The AR(2)-ARCH(4) model was certainly an improvement on ARMA(2,1), judging from the values
of model selection criteria. However, as we know from the LM test result, including only four lags
of ✏2
t in the specification of ht seems inadequate. Including more lags in (1) and re-testing reveals
that ˆ ✏2
t from ARMA(2,1) model in fact depend on ˆ ✏2 t from the past two weeks or more. In your
own time, conduct the LM test for ARCH errors with varying lag lengths. Here is an example with
12 lags (recall that this is daily data, so 12 lags is just under two-and-a-half weeks of trading).
7
It has been reported in the literature that one often needs to include a very long lag in ARCH
models to capture the changing volatility adequately, leading to many parameters to estimate. It
has also been noted that a GARCH model of very small orders fits the data just as well as an
ARCH model of a high order. We shall now estimate a GARCH(1,1) model on the same data, and
see if it gives us a better fit to the data. We shall discard the two AR components as they both
turned out to be totally insignificant in the above AR(2)-ARCH(4) estimation. The model we will
fit is specified as follows:
yt = c + ✏t (6)
where

↵0 > 0, ↵1 ” 0, &1 ” 0 and ↵1 + &1 < 1
In the main gretl window, click on Model #! Time series #! GARCH… again to bring up the
GARCH estimation dialogue box. Select r tm for the dependent variable, no regressors, and select
1 for both GARCH p: and ARCH q: as shown on the next page.
8
Click on OK to obtain the following estimation output.
We observe further improvements in all three criteria, although not to the same extent as we
observed in AR(2)-ARCH(4) model from ARMA(2,1). Notice that &1 is suspiciously significant,
and that ↵1 + &1 ⇡ 1. This is very often observed when GARCH(m, r) model is fitted to real
financial data—the requirement for stationarity of ht, which is
mX i
=1
↵i +
rX j
=1
&j < 1 (9)
is very close to being violated. When condition (9) is violated—in the simplest case of GARCH(1,1),
↵1 + &1 = 1—such a process is called IGARCH, for integrated GARCH. Whether or not such a
process makes sense for an obviously stationary financial return series is a matter of debate.
9
In your own time, take a look at the estimated ht from this model. Compare it with the one
obtained from AR(2)-ARCH(4) above. Also, create the standardised residuals from GARCH(1,1):
v˜tg =
✏˜tg
qh˜tg
where the ‘g’ in subscript denotes that the residuals, the estimated ht, and the standardised residual
are from the GARCH(1,1) model. Plot the correlogram of ˜ vtg. Here is what you should obtain:
Does the above graph indicate an improvement from fitting AR(2)-ARCH(4) to r tm? Create the
squared standardised residuals ˜ vtg 2 and plot its correlogram.
3 Simulating GARCH processes in gretl
Ensure that all the estimation and testing outputs you have so far have been saved as icons, then
save any changes to the session file by clicking on File, then Save data, as well as on Save session.
We will be using the nulldata command again.
As another simulation exercise, we will simulate a GARCH(1,1) process as defined in (6), (7)
and (8). We shall first simulate the process using the parameter estimates obtained from fitting a
GARCH(1,1) model on Microsoft stock returns:
c = 0.000885681
↵0 = 0.00000524367
↵1 = 0.062619
&1 = 0.928127
h1 = 0.00056662
where h1 denotes the initial condition. In order to start a loop to create (8), the value for h1,
the conditional volatility for the very first time period, needs to be specified. We shall use the
unconditional error variance estimated on Microsoft stock returns for this.
10
In any simulations, all of which typically depend on initial conditions, it is wise to discard
an early proportion of the series which could well su↵er from the e↵ect of the initial conditions.
Therefore, we shall simulate a series of length 5,000, but discard the first 1,000 which will give us
a series of length 4,000.
Open the script file GARCH 11 sim.inp that I gave you. It has declarations of the parameters,
the sample range, as well as lengthy comments on what the script does. Starting from line 34, type
the following:
34 series e t=zeros(n,1)
35 series h t=zeros(n,1)
36 v t=mnormal(n,1)
37 h t[1]=h 1
38 e t[1]=c+sqrt(h 1)*v t[1]
39
40 loop i=2..n –quiet
41 h t[i]=alpha 0+alpha 1*(e t[i-1])^2+beta 1*h t[i-1]
42 e t[i]=c+sqrt(h t[i])*v t[i]
43 endloop
44 smpl 1001 n
45
46 p t=cum(e t)
47 series Price=exp(p t)
Some of the commands should be familiar to you from the last class. The explanations for the
above lines are given below.
34. and 35. Initialise two n⇥1 null vectors to store simulated values of ✏t and ht. e t corresponds
to the observed series ✏t (plus constant, see point 42 below) and h t is to store the simulated
conditional volatility series ht.
Note that they are both preceded by command series. This declares e t and h t to be
vector series, and not matrices. This is essential when you wish to plot a simulated series.
Even if you simulate a 10 ⇥ 1 series which is obviously a vector, unless it is preceded by the
command series, it is a matrix as far as gretl is concerned. When gretl recognises a time
series as a matrix, it will be much more complicated to plot it. Whatever you declare as a
series, gretl automatically stores it as data in the main gretl window. Needless to say, every
series needs to be of the length declared after the command nulldata.
36. v t=mnormal(n,1): Takes n independent draws from the standard normal distribution and
stores it in an n⇥1 vector. This is vt in equation (7). Note that v t is a matrix as far as gretl
is concerned as we have not declared it a series. Of course you can precede this command
by series—then gretl will recognise it as a vector of series, and will store it as part of the
dataset. But since this is not a time series of interest which we wish to plot, I wouldn’t
bother.
Also note that whenever you store or print anything at all, this adds to the computational
time as well as taking up storage space. In this instance, we are simulating one time series,
so the e↵ect of any additional time or space is negligible. However, in a simulation exercise
involving simulating a million series, it will make a considerable di↵erence.
11
37. h t[1]=h 1: Stores the value of h 1, declared on line 32 as a scalar, into the first element of
h t. Here, we are specifying the starting value of the loop, h1.
38. e t[1]=c+sqrt(h 1)*v t[1]: Once h1 is known, this line of command specifies ✏1, in e t[1],
according to equation (6), and adds the constant. Here, command e t[1]=c+sqrt(h 1)*v t[1,1]
will also be acceptable to gretl, because v t is a matrix as far as gretl is concerned. A matrix
has two dimensions, so a specific element within it is specified by its row and column locations. However, the command e t[1,1]=c+sqrt(h 1)*v t[1] will be an incorrect syntax
and gretl will fail to execute it—this is because e t has been declared as a vector of series.
40. loop i=2..n –quiet: Declares the beginning of a loop, telling gretl that the counter i
needs to begin from 2 and ends at n. Note that the first elements of e t and h t have already
been filled. Therefore we wish to fill the remainder of them, generating values according to
(6), (7), and (8). –quiet supresses printing thereby speeding up the computation.
41. Fills h t as specified in (8). This line specifies ht = ↵0 + ↵1✏2 t!1 + &1ht!1.
42. Fills e t as specified in (7), then adds the constant to create (6).
Note two things here: one, in every round of the loop, that is for a particular i, the loop
fills the ith element of both e t and h t. It does not (and cannot) fill up all elements of one
vector before starting to fill the other. Two, the equation for ht needs to be specified before
that of ✏t. This is because for every i, to fill the ith row of e t, the ith row of h t must
already be filled, as ✏t = phtvt.
43. endloop: Self-explanatory.
44. smpl 1001 n: Resets the sample range. Here we are discarding the first 1,000 simulated
values, and only saving the last 4,000.
46. p t=cum(e t): Takes the cumulative sum of e t. Recall that GARCH(1,1) was fitted on the
log returns on Microsoft stocks. This was created by taking the first di↵erence of the natural
log of the raw data on Microsoft stock price. Taking the cumulative sum of the di↵erenced
series will give us back the original level series. So p t corresponds to the log of stock price.
47. series Price=exp(p t): Takes the exponent of every element of p t, so essentially unlogging the logged series p t. So Price corresponds to the raw data on stock price, and is
declared as a series, enabling us to plot it easily.
Once you have completed typing in all the commands, save it by clicking on the floppy disc icon
in the toolbar. Execute the code by clicking on the cogwheel icon.
You should see several data series in your gretl window, including an index. These are the
variables that have been declared as series or have been created from declared series. Click on
the icon for gretl console, the third from the left at the bottom of the main gretl window. At the
prompt, type the following command:
gnuplot Price –time-series –with-lines
and hit the return key. The above command plots series Price as time series using a line graph.
Every time you plot any graph, whether it be time series or histogram, gretl calls another opensource software package called gnuplot. This is a stand-alone software which was developed to
12
visualise data1. The default for plotting a series in gnuplot is using a cross marker to mark the
location of every observation. Since line graph is much easier to see, –with-lines is an argument
permitted in command gnuplot to change the default to a line graph. You should obtain the
picture below, plotted on the left (provided that you have not changed the seed of the random
number generator from 12345).
Simulated Price and ht with &1 = 0.928127
Back in the window for gretl console, press the upward arrow key on your keyboard. It should
bring up the last command you typed on the next commmand line. Now we wish to visualise the
simulated ht series. So change the Price in the line of command you have just brought up to h t
to produce the command below:
gnuplot h t –time-series –with-lines
and hit the return key. You should obtain the graph on the right-hand-side above. Compare it
with the fitted ht from estimation of GARCH(1,1) on Microsoft stock returns.
We observed earlier that the GARCH(1,1) model fitted on Microsoft returns—which we are
simulating—is very close to being an integrated GARCH. In the script file you have just completed,
change the value of beta 1 from 0.928127 to 0.5 (Note: it’s much easier to comment out line 31
by inserting a hashtag in front of it, then writing scalar beta 1=0.5 on the next line. This way,
when you wish to restore the original &1, you don’t have to type it in.) Execute the script.
Now we wish to plot Price and h t. Back on the gretl console, press the upward arrow key on
your keyboard twice. It should bring up the line of command you typed in not the last time, but the
one before that, which was gnuplot Price –time-series –with-lines. This is exactly what
we want, so hit the return key. You should obtain a much smoother upward sloping curve plotted
on the next page. In gretl console, press the upward arrow twice again to bring up gnuplot h t
–time-series –with-lines and hit the return key. The picture you should obtain is plotted
on the right-hand-side below. Compare to above where &1 = 0.928, the conditional volatility is far
more even across time.
In GARCH 11 sim.inp script, leave &1 at 0.5 and now set constant c as 0. Execute the code
again. Using gretl console as before, plot series Price and h t again. They are printed on the next
page as well.
1gnuplot can be downloaded from here: http://gnuplot.sourceforge.net/, and like gretl, can be installed legally
for free on any computer
13
Simulated Price and ht with &1 = 0.5
Simulated Price and ht with &1 = 0.5 and c = 0
We can observe the role the constant term was playing, although it was miniscure at 0.000886.
Since this was part of an exponent in Price, the cumulative e↵ect was significant on the upward
exponenrial trend of the process. When &1 was close to 1, the sheer change in volatility of the series
seems to have had a dampening e↵ect on the exponential growth of Price. Note that the series
h t is nearly identical to the one plotted with a positive constant, for an obvious reason that the
parameters governing ht have not been altered this time.
Experiment with di↵erent values for all the parameters of the process, and see what di↵erence
it makes to the resulting series. Also plot simulated e t and p t.
4 Estimating non-linear GARCH variants
Estimate a variety of GARCH variant models in your own time, and see if there is any evidence of
leverage e↵ect and asymmetry in series r tm.
As discussed in the lectures, one weakness of both ARCH and GARCH models is symmetry—
only the magnitude of past innovations influence the current conditional volatility, and the sign of
the past shocks have not been taken into consideration. Therefore, positive and negative shocks of
the same magnitude are assumed to have the same influence on current ht. However, it has often
been observed that a negative shock a↵ects the volatility of financial time series more profoundly
than a positive shock of the same magnitude. Exponential GARCH developed by Nelson (1991)
is one of the models which permit asymmetry between the impact of positive and negative shocks
14
on conditional volatility. It also circumvents the non-negativity conditions required for ARCH and
GARCH models by modelling the natural log of ht.
In its simplest form, EGARCH(1,1) is defined as follows:
✏t = c + phtvt; where vt ⇠ iid(0, 1) (10)
where
ln(ht) = ↵0 + g(vt!1) + &1 ln(ht!1)
g(vt) = ✓vt + ( [|vt| # E(|vt|)]
By substituting g(vt) into the equation for ln(ht) and rearranging, the conditional variance equation
is given by
ln(ht) = ↵0⇤ + ✓vt!1 + (|vt!1| + &1 ln(ht!1) (11)
where ↵⇤
0 = ↵0#(p2/⇡ since E(|vt!1|) = E(|vt|) = p2/⇡ ⇡ 0.8. This is the equation estimated by
gretl. However, they use di↵erent notations for the parameters. In their notation, EGARCH(1,1)
gretl estimates is defined as:
ln(ht) = ! + ↵|vt!1| + (vt!1 + & ln(ht!1) (12)
so gretl’s ↵ is ( in my (and Nelson’s original) notation, its ( is our ✓. This is the parameter which
captures the asymmetry in the impact positive and negative shocks have on volatility. If, as often
observed, negative shocks produce greater volatility while positive shocks dampen the series, we
expect ( (our ✓) to be negative.
gretl has a package called gig built in to estimate a variety of GARCH-type models. Clicking on Models #! Time series #! GARCH variant… gives you the dialogue box shown on the
left. There is a drop-down menu where you can
choose which GARCH variant you wish to estimate. As you can see, the most common models
such as threshold GARCH, GJR-GARCH, and
EGARCH, are available amongst others. You
can select a model, the number of lags, click on
OK to estimate a model. However, using this
menu-driven way will not give you an option to
save the residuals or the estimated ht. In order
to do so, you need to call gig manually through
gretl’s console or by containing commands in a
script.
The script file egarch.inp is provided on our
Blackboard site for this purpose. I have written
a few lines to call gig, save the residual series as
e and the estimated ht as h. Open egarch.inp,
and execute it by clicking on the cogwheel icon.
You should obtain the following in your script
output (estimation of this model may take several seconds).
15
The estimated ( is negative as expected, and statistically significant at 5% level, given a p-value
of just under 0.02. Therefore it appears that negative shocks do have a greater impact on ht than
positive ones. The coe”cient of ln(ht!1) is very close to 1, with a suspiciously high t-ratio—this
is consistent with the estimation result of GARCH(1,1) in which the estimated coe”cients were
bordering on the limits of the stationarity condition for ht.
Also notice that the values for all three criteria have been reduced yet further from those we
obtained from GARCH(1,1). This indicates that taking asymmetry into consideration does indeed
improve the fit of the model.
You can save the entire content of the script output by clicking on the floppy disc icon in the
script output window, and selecting save to session as an icon. You should see script output saved in
your session icon view. You should also see a folder named emod. This contains the EGARCH(1,1)
estimation result. Note that in egarch.inp script, I had written a line of code to print the sample
variance of fitted ht and its mean. These will be saved in the script output, but not in the estimation
output saved in emod folder.
In your main gretl window, you should see two new series added: h and e. They are the fitted
conditional volatility series obtained from (12), and the residual series ˜ ✏t respectively. h is given
by:
h = exp hln(h˜t)i (13)
where ln(h˜t) is given by
ln(h˜t) = ˜ ! + ˜ ↵|v˜t!1| + ˜ (v˜t!1 + &˜ln(h˜t!1)
where ˜ denotes the ML estimation of parameters, and fitted values using these estimates. Using
the estimates of conditional volatility and the residuals, construct the standardised residuals, ˜ vt.
Plot its correlogram. You should obtain the following:
16
Does the above correlogram indicate a further improvement in the fit of the model from GARCH(1,1)?
Take a look at the correlogram of ˜ vt2 as well.
5 Further exercises
1. Estimate a variety of ARCH, GARCH, and other variants of GARCH models on the return
series of Apple and HSBC stocks contained in your dataset. See if you observe any evidence
for asymmetry and leverage e↵ect.
2. In the GARCH simulation script, change the parameters of the process. Simulate GARCH
process ✏t and conditional volatility series ht for a variety of parameter combinations. What
relationship do you observe between the values of any one of the parameters (and all of the
parameters as a set) and the behaviour of the resulting ✏t and ht? When the stationarity
condition ↵1 + &1 is violated, what does the series look like?
3. Try to simulate EGARCH(1,1) process. Try a variety of values for ( (that’s gretl’s (), the
parameter which captures asymmetry.
17