```----------------------------------------------------
#                                                  #
#    Chapter 5: Maximum Likelihood Estimation      #
#                of the ARCH(1)-Model              #
#                                                  #
----------------------------------------------------

In the last chapter we defined the ARCH(d)-model and
we simulated a couple of paths for various choices
of model paramters. We saw that, depending on the
choice of these model parameters, the model is already
able to generate quite a variety of different price
dynamics.

Now the question arises how we can determine which of
these price dynamics, that is, which choice of model
parameters, are the best for a given time series.
In this chapter we will answer this question for the
ARCH(1)-model. The logic for the ARCH(d)-model is
completely analog and will be considered in the next
chapter.

The ARCH(1)-Model is given by

S(t_k) = S(t_{k-1}) * [ 1 + vol(t_{k-1})*phi_k ]                  (1)

with volatility specification

vol^2(t_{k-1}) = w0*bsvol^2 + w1*ret^2(t_{k-1})                   (2)

with weights  w0+w1=1  and a constant Black-Scholes
volatility bsvol. For  w0=1, w1=0  the model reduces
to the time discrete Black-Scholes model (with zero
drift).

Let's fit this model to the DAX time series. That is,
we would like to determine the parameters w0 (then
w1=1-w0 is also fixed) and bsvol, such that the
corresponding price dynamics fits `best' to the given
time series data. The meaning of `best' of course
has to be formulated in some precise mathematical way.

To this end we rewrite (1) as

ret(t_k) = vol(t_{k-1}) * phi_k                                   (3)

and recall, that phi_k is a standard normal-distributed
random number. Thus, if we are at time t_{k-1}, then
the probability that the next return ret(t_k) will
realize in the intervall [y_k,y_k+dy_k], is given by

Prob[ ret(t_k) in [y_k,y_k+dy_k] | t_{k-1} ]                      (4)

= 1/sqrt{2pi} 1/vol(t_{k-1}) * exp( -1/2*y_k^2/vol^2(t_{k-1}) ) dy_k

Thus, the probability that all returns have values in
some intervalls [y_k,y_k+dy_k], for all k=1,2,...,N,
is given by the product of the probabilities in (4),

Prob[ DAX-ret in [y_k,y_k+dy_k] for all k ]

= prod_{k=1}^N 1/sqrt{2pi} 1/vol(t_{k-1}) *                     (5)

prod_{k=1}^N exp( -1/2*y_k^2/vol^2(t_{k-1}) ) dy_k

Now, if we substitute the ARCH(1) vol-specification given
by (2) above,

vol^2(t_{k-1}) = w0*bsvol^2 + w1*y_{k-1}^2                        (2)

on the right hand side of (5), then the probability in (5),

Prob[ DAX-ret in [y_k,y_k+dy_k] for all k ]

can be considered as a function of the model parameters
bsvol and w0 exclusively, since the y_k are actually to be
substituted by the actual realized returns of the last 10
years. To make this more explicit, that we actually consider
this probability as a function of the model parameters
bsvol and w0, we introduce a new letter L (for "Likelihood")
and write

Prob[ DAX-ret in [y_k,y_k+dy_k] for all k ]

=  L( w0 , bsvol )                                              (6)

=  L( w0 , bsvol | y_k=ret(t_k) )

where the notation in the last line of (6) should just be a
reminder that the y_k are given data, namely the returns
which have realized. This L is then called the likelihood
function and a standard procedure in statistics is simply to
estimate model parameters by maximizing this likelihood func-
tion. Accordingly, this method is called "Maximum Likelihood
Estimation".

Since we have about N=2500 returns in a time series of 10
years, L is given by a product of N=2500 factors. Thus, we
should expect that L behaves like

L ~ exp( +N*something )  or   L ~ exp( -N*something )           (7)

which is either extremely large or extremely small. To avoid
numerical issues, it is therefore common to consider the
logarithm log(L) of the likelihood function. Apparently,
maximizing L is equivalent of maximizing log(L), so the logic
stays the same, but for log(L) we can expect a behaviour of

log(L) ~  +N*something   or   log(L) ~  -N*something            (8)

which is much more tractable than the expressions in (7).
Thus we consider

log(L) = sum_{k=1}^N { log[1/vol(t_{k-1})] -                      (9)

1/2 * y_k^2 / vol^2(t_{k-1}) }
+ const

where the constant

const = N*log[1/sqrt{2pi}] + sum_{k=1}^N log[dy_k]               (10)

does not depend on the model paramters w0 and bsvol. Since
we maximize with respect to model parameters, we therefore
can ignore the constant (10) and end up with the problem of
maximizing the function

log(tilde L) = log(tilde L)( bsvol , w0 )                        (11)

= - sum_{k=1}^N { log[vol(t_{k-1})] + 1/2 * y_k^2 / vol^2(t_{k-1}) }

with  y_k = ret(t_k)  and

vol^2(t_{k-1}) = w0*bsvol^2 + (1-w0)*ret^2(t_{k-1})               (2)

This maximization problem has no longer an analytic closed
form solution and one has to resort to some kind of numerical
method. This is actually a typical situation in Maximum Likeli-
hood Estimation: the maximum cannot be found in analytic closed
form, but has to be computed numerically.

We make an R-implementation to perform this task:

# Start R-Session:
#
# before we fit to DAX-data, we check on an artificially
# generated ARCH(1) return set: let's use the pathArchd
# function from the last chapter with d=1 to generate these
# artificial return data:

# to make this reproducible, we fix a seed for
# the random number generator:
set.seed(456)

pathArchd = function( N , bsvol , w0 , d )
{
S = rep(0,N)
ret = rep(0,N)
phi = rnorm(N)
S0 = 100
vol0 = bsvol            # we put startvol to bsvol
vol = rep(vol0,N)       # ARCH-vol

sumretd = 0

# the first step, k=1, is special:
S[1] = S0*(1+vol0*phi[1])
ret[1] = (S[1]-S0)/S0
vol[1] = sqrt( w0*bsvol^2 + (1-w0)*ret[1]^2 )
sumretd = sumretd + ret[1]^2

for(k in 2:N)
{
if(k<=d)
{
S[k] = S[k-1]*(1+vol[k-1]*phi[k])
ret[k] = (S[k]-S[k-1])/S[k-1]
sumretd = sumretd + ret[k]^2
vol[k] = sqrt( w0*bsvol^2 + (1-w0)*sumretd/k )
}
else
{
S[k] = S[k-1]*(1+vol[k-1]*phi[k])
ret[k] = (S[k]-S[k-1])/S[k-1]
sumretd = sumretd + ret[k]^2
sumretd = sumretd - ret[k-d]^2
vol[k] = sqrt( w0*bsvol^2 + (1-w0)*abs(sumretd)/d )
}
}#next k

par(mfrow=c(2,1))
plot(S,type="l",ylim=c(0,250))
plot(ret,type="l",ylim=c(-0.08,0.08))

result = list(S,ret)
names(result) = c("S","ret")
return(result)
}

N = 2500
bsvol = 0.01
w0 = 0.5
res = pathArchd( N , bsvol , w0 , d=1 )
ret=res\$ret
summary(ret)

# artificially generated price path
# and return set:
```

```
# for this return set ret we now want to calculate the
# log-likelihood function and a maximization should
# result in model parameters around bsvol = 1% = 0.01
# and w0 = 0.5:
#
# we calculate the log-Likelihood function:

logL = function( bsvol , w0 )
{
w1 = 1 - w0
vol = sqrt(w0*bsvol^2+w1*ret^2)     # we can use the vector ret, although we
# did not pass it through the arguments
# shift vol by 1 element:
n = length(ret)
vol = vol[-n]                       # remove last one
vol = c(sqrt(w0*bsvol^2),vol)       # add a new first one
terms = log(vol) + 1/2*ret^2/vol^2
result = -sum(terms)
return(result)
}

# let's check:
logL( bsvol=0.01 , w0=0.5 )           # ok, we get some number

# since we have only 2 variables, we can look for the maximum
# by making a contour-plot. To become familiar with the
# contour()-function, let's first consider an example:

#-----------------------------------------------------
# START Example: Contour-Plot:
#
# function of 2 variables to be plotted:

testf = function(x,y)
{
result = (-x^2+y^2)*exp(-0.25*x^2-0.5*y^2)
return(result)
}

# the contour()-command has syntax
# contour(  vector ,  vector , matrix ) =
# contour( x-Achse , y-Achse , f(x,y) )

x = seq(from=-4,to=4,by=0.1)
y = x

# matrix with testf(x,y)-values:
z = outer( x , y , FUN = testf )     # this works only if the function testf can
# accept vectors as input variables

# and now the level curves:
contour( x , y , z )
```

```
contour( x , y , z , nlevels=50 )    # more refined
```

```
# instead of using the outer()-command
# we can use a loop to generate z:

nx = length(x)
ny = length(y)
z = matrix( 0 , nrow=nx , ncol=ny )
for(i in 1:nx)
{
for(j in 1:ny)
{
z[i,j] = testf(x[i],y[j])        # now scalars as arguments,
# no longer vectors
}
}

# same contour-plots as above:
contour( x , y , z )
contour( x , y , z , nlevels=50 )    # more refined

# also useful:
image( x , y , z )
```

```
persp( x , y , z )
```

```
# with plot3D package:
require(plot3D)
persp3D( x , y , z )
```

```#
# END Example Contour-Plot
#-----------------------------------------------------

# Let's make a contour-plot of the log-likelihood function:
#
# Grid for the model parameters (recall that the model
# is unstable for w0=0 or bsvol=0):

bsvol = seq( from=0.1/100 , to=4/100 , by=0.1/100 )
w0 = seq( from=0.05 , to=0.95 , by=0.05 )

# the following does not work, since logL does not accept
# vectors as input; the reason is that in logL there is a line
# vol = sqrt(w0*bsvol^2+w1*ret^2)
# which does not make sense if all w0, w1 and ret are vectors:

z = outer( bsvol , w0 , FUN = logL )          # does not work

# so let's do this with a loop:

nbs = length(bsvol)
nw = length(w0)
z = matrix( 0 , nrow=nbs , ncol=nw )
for(i in 1:nbs)
{
for(j in 1:nw )
{
z[i,j] = logL( bsvol[i] , w0[j] )         # this is ok
}
}

# let's look at the result:
contour( bsvol , w0 , z )
```

```
contour( bsvol , w0 , z , nlevels = 50 )
```

```
contour( bsvol , w0 , z , nlevels = 50 , zlim=c(10000,11000) )
# eventually vary range for zlim-values
```

```
# okay, that works!!  -> we recover correct parameter values!!
# bsvol around 0.01 and w0 around 0.5

# some more pictures (would require more fine-tuning):
image( bsvol , w0 , z )
```

```
image( bsvol , w0 , z , col = rainbow(100) , zlim=c(10000,11000) )
```

```
persp( bsvol , w0 , z )
```

```
persp3D( bsvol , w0 , z )
```

```

---------------------------------------
#     ARCH(1)-Fit to DAX-Data:        #
---------------------------------------

# load DAX.txt data and generate return vector ret:

dax = na.omit(dax)
summary(dax)
S = dax\$index
n = length(S)
ret = rep(0,n)
for(i in 2:n)
{
ret[i] = (S[i]-S[i-1])/S[i-1]
}
plot(ret,type="l",main="dax-returns")
```

```
#
# calculate DAX-logL and make setup for contour-plot:
#
bsvol = seq( from=0.1/100 , to=4/100 , by=0.1/100 )
w0 = seq( from=0.05 , to=0.95 , by=0.05 )
nbs = length(bsvol)
nw = length(w0)

z = matrix( 0 , nrow=nbs , ncol=nw )
for(i in 1:nbs)
{
for(j in 1:nw )
{
z[i,j] = logL( bsvol[i] , w0[j] )
}
}

# let's look at the result:
contour( bsvol , w0 , z )
```

```
contour( bsvol , w0 , z , nlevels = 50 )
```

```
contour( bsvol , w0 , z , nlevels = 50 , zlim=c(8000,10000) )
```

```
contour( bsvol , w0 , z , nlevels = 50 , zlim=c(9000,9600) )
```

```# we take a more refined look at the region w0 = 0.5 to 0.9
# and bsvol = 1.0% to 2.0%:

bsvol = seq( from=1.0/100 , to=2.0/100 , by=0.02/100 )
w0 = seq( from=0.5 , to=0.9 , by=0.01 )
nbs = length(bsvol)
nw = length(w0)
z = matrix( 0 , nrow=nbs , ncol=nw )
for(i in 1:nbs)
{
for(j in 1:nw )
{
z[i,j] = logL( bsvol[i] , w0[j] )
}
}

# let's look at the results:
contour( bsvol , w0 , z , nlevels = 100 , zlim=c(9000,9600) )
```

```contour( bsvol , w0 , z , nlevels = 100 , zlim=c(9500,9550) )
```

```max(z)                             # 9534.98
logL( bsvol=sd(ret) , w0=0.72 )    # 9534.64
sd(ret)                            # 1.39%

# thus: maximum approximately at
# bsvol = sd(ret) = 1.39% and w0 = 0.72

# Before we redo the calculations for SPX and GE, we code
# a small function GetZ which sets up the matrix z for
# the contour-plots:

GetZ = function( x , y )
{
nx = length(x)
ny = length(y)
z = matrix( 0 , nrow=nx , ncol=ny )
for(i in 1:nx)
{
for(j in 1:ny)
{
z[i,j] = logL(x[i],y[j])
}
}
return(z)
}

---------------------------------------
#     ARCH(1)-Fit to S&P500-Data:     #
---------------------------------------

# load SPX.txt data and generate return vector ret:

summary(spx)     # no NA's
S = spx[,2]
n = length(S)
ret = rep(0,n)
for(i in 2:n)
{
ret[i] = (S[i]-S[i-1])/S[i-1]
}
plot(ret,type="l",main="spx-returns")   # returns look ok
```

```
# compute spx-logL and make setup for contour-plot:

bsvol = seq( from=0.1/100 , to=4/100 , by=0.1/100 )
w0 = seq( from=0.05 , to=0.95 , by=0.05 )

z = GetZ( bsvol , w0 )

# let's look at the results:

contour( bsvol , w0 , z )
```

```
contour( bsvol , w0 , z , nlevels = 50 )
```

```contour( bsvol , w0 , z , nlevels = 50 , zlim=c(60000,80000) )
```

```contour( bsvol , w0 , z , nlevels = 100 , zlim=c(60000,80000) )
```

```# we take a more refined look at the region w0 = 0.5 to 0.9
# and bsvol = 0.5% to 1.5%:

bsvol = seq( from=0.5/100 , to=1.5/100 , by=0.02/100 )
w0 = seq( from=0.5 , to=0.9 , by=0.01 )

z = GetZ( bsvol , w0 )

contour( bsvol , w0 , z , nlevels = 100 , zlim=c(60000,80000) )
```

```contour( bsvol , w0 , z , nlevels = 100 , zlim=c(69000,70000) )
```

```max(z)                             # 69471.9
logL( bsvol=sd(ret) , w0=0.69 )    # 69472.17
sd(ret)                            # 0.97%

# thus: maximum approximately at
# bsvol = sd(ret) = 0.97% and w0 = 0.69

---------------------------------------
#     ARCH(1)-Fit to GE-Data:         #
---------------------------------------

# load GE.txt data and generate return vector ret:

summary(ge)     # no NA's
S = ge[,2]
n = length(S)
ret = rep(0,n)
for(i in 2:n)
{
ret[i] = (S[i]-S[i-1])/S[i-1]
}
plot(ret,type="l",main="General Electric returns")
# returns look ok:
```

```# compute GE-logL and make setup for contour-plot:

bsvol = seq( from=0.1/100 , to=4/100 , by=0.1/100 )
w0 = seq( from=0.05 , to=0.95 , by=0.05 )

z = GetZ( bsvol , w0 )

# let's look at the result:

contour( bsvol , w0 , z )
```

```contour( bsvol , w0 , z , nlevels = 50 )
```

```contour( bsvol , w0 , z , nlevels = 50 , zlim=c(40000,50000) )
```

```# we take a more refined look at the region w0 = 0.5 to 0.9
# and bsvol = 1.0% to 2.0%:

bsvol = seq( from=1.0/100 , to=2.0/100 , by=0.02/100 )
w0 = seq( from=0.5 , to=0.9 , by=0.01 )

z = GetZ( bsvol , w0 )

contour( bsvol , w0 , z , nlevels = 50 , zlim=c(40000,50000) )
```

```contour( bsvol , w0 , z , nlevels = 100 , zlim=c(49500,50000) )
```

```max(z)                             # 49811.81
logL( bsvol=sd(ret) , w0=0.68 )    # 49811.89
sd(ret)                            # 1.62%

# thus: maximum approximately at
# bsvol = sd(ret) = 1.62% and w0 = 0.68

```