```---------------------------------------------------
#                                                 #
#     Chapter 6: Maximum Likelihood Estimation    #
#              of the ARCH(d)-Model               #
#                                                 #
---------------------------------------------------

In the last chapter we familiarized ourselves with the
maximum likelihood method and we derived the following
expression for the log-likelihood function (we discarded
a constant involving logarithms of 2*pi's):

log(L) =                                                          (1)

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

The y_k are given by the realized returns of the actual
time series data,

y_k = ret(t_k) = [ S(t_k)-S(t_{k-1}) ] / S(t_{k-1})               (2)

and the exact form of the volatility vol(t_k) depends
on the model under consideration. For the ARCH(1)-model,
it is given by

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

with w0 and bsvol being the model parameters. For the
DAX time series, we found

DAX:    w0 = 0.7  und  bsvol = 1.4%                               (4)

by maximizing the ARCH(1) log-likelihood function. The
parameter bsvol of 1.4% turned out to be very close to
the standard deviation sd(ret) of the returns taken over
the whole time series length, sd(ret) = 1.39%.

The more realistic ARCH(d)-model which is to be estimated
in this chapter has the following vol-specification:

vol^2(t_{k-1}) =

w0*bsvol^2 + (1-w0)*[ret^2(t_{k-1})+...+ret^2(t_{k-d})] / d      (5)

with ARCH(d) model parameters w0, bsvol and the time
horizon d. The maximization of the likelihood function
again has to be done numerically and to this end we
start an R-session:

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

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)
}

# we generate an artificial ARCH(d) return set where we
# know the actual model parameters ( bsvol , w0 , d ):
# let's set a seed for the random numbers to make this
# reproducible:
set.seed(234)

N = 2500
bsvol = 0.01
w0 = 0.3
d = 20
res = pathArchd( N , bsvol , w0 , d )
```

```ret = res\$ret
summary(ret)

# we have to calculate the log-Likelihood function:
#
# use 'vectorized' calculation logic of R to make the code
# a bit more performant -> code looks a bit less intuitive:

logL = function( bsvol , w0 , d )
{
w1 = 1 - w0
n = length(ret)

cumsum_retsquared = cumsum(ret^2)

# remove last d elements:
cumsum_retsquared_shifted = cumsum_retsquared[-((n-d+1):n)]

# fill in d zeros at the first d places:
cumsum_retsquared_shifted = c(rep(0,d),cumsum_retsquared_shifted)

sum_dretsquared = cumsum_retsquared - cumsum_retsquared_shifted

# now we have to divide this sum by d, but the first d elements
# are special:
divisor1 = 1:d
divisor2 = rep(d,n-d)
divisor = c(divisor1,divisor2)

sum_dretsquared = sum_dretsquared/divisor

# now we can calculate ARCH(d)-vol:
vol = sqrt( w0*bsvol^2 + w1*sum_dretsquared )

# shift vol-vector by 1 element ( there is a vol(t_{k-1}) in
# the likelihood-function, not a vol(t_k) )
vol = vol[-n]                           # remove last one
vol = c(sqrt(w0*bsvol^2),vol)           # add a new first one

# the actual calculation:
terms = log(vol) + 1/2*ret^2/vol^2
result = -sum(terms)

return(result)
}

# let's check:
logL( bsvol=0.01, w0=0.3, d=20 )          # alright, some number

# let's look for the maximum:
# first, we fix bsvol to its natural value:
bsvol = sd(ret)
bsvol

# we make a contour-plot in the (w0,d)-plane:

w0 = seq(from=0.05,to=0.95,by=0.05)
dd = 1:40
nw = length(w0)
nd = length(dd)

# matrix z = logL for contour-plot:
z = matrix( 0 , nrow=nd , ncol=nw )
for(i in 1:nd)
{
for(j in 1:nw )
{
z[i,j] = logL( bsvol , w0[j] , dd[i] )
}
}# still very fast evaluation of logL!

# let's look at the result:
contour( dd , w0 , z , nlevels = 50 )
```

```contour( dd , w0 , z , nlevels = 50 , zlim=c(10000,11000) )
```

```contour( dd , w0 , z , nlevels = 100 , zlim=c(10200,10600) )
```

```
# we recover correct parameter values:
# d around 20 and w0 around 0.3;
# differences between d = 15, 20 and 25 not that big,
# depending on the particular realization of random numbers;
# d below 8 clearly less favorable -> redo calculation for
# dd = 10:40 and w0 < 0.5 (-> resolution on pictures more nice)

w0 = seq(from=0.05,to=0.5,by=0.01)
dd = 10:40
nw = length(w0)
nd = length(dd)
z = matrix( 0 , nrow=nd , ncol=nw )
for(i in 1:nd)
{
for(j in 1:nw )
{
z[i,j] = logL( bsvol , w0[j] , dd[i] )
}
}
contour( dd , w0 , z , nlevels = 150 , zlim=c(10200,10600) )
```

```# some more pictures:
image( dd , w0 , z )
```

```persp( dd , w0 , z )
```

```require(plot3D)
persp3D( dd , w0 , z )
```

```
# let's fix d=20 and check that bsvol=sd(ret) corresponds
# to the maximum:

bsvol = seq(from=0.1/100,to=4/100,by=0.1/100)
nbs = length(bsvol)
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] , d=20 )
}
}
contour( bsvol , w0 , z , nlevels = 50 )
```

```contour( bsvol , w0 , z , nlevels = 80 , zlim=c(10000,10500) )
```

```# -> looks all good.

#
# finally, let's fix w0 and vary bsvol and d:
#
bsvol = seq(from=0.1/100,to=2/100,by=0.05/100)
nbs = length(bsvol)
dd = 10:40
nd = length(dd)
z = matrix( 0 , nrow=nd , ncol=nbs )
for(i in 1:nd)
{
for(j in 1:nbs )
{
z[i,j]=logL( bsvol[j] , w0=0.3 , d=dd[i] )
}
}
contour( dd , bsvol , z , nlevels = 50 )
```

```

-----------------------------------------
#      ARCH(d)-Fit to DAX-Data:         #
-----------------------------------------
dax = na.omit(dax)
summary(dax)
rm(ret)
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")
```

```
w0 = seq(from=0.05,to=0.95,by=0.05)
dd = 1:40
bsvol = sd(ret)
bsvol
nw = length(w0)
nd = length(dd)
z = matrix( 0 , nrow=nd , ncol=nw )
for(i in 1:nd)
{
for(j in 1:nw )
{
z[i,j] = logL( bsvol , w0[j] , dd[i] )
}
}

# let's look at the result:
contour( dd , w0 , z , nlevels = 50 )
```

```contour( dd , w0 , z , nlevels = 50 , zlim=c(9500,10500) )
```

```# apparently w0 less than 0.4 and d between 8 and 28:

w0 = seq(from=0.01,to=0.4,by=0.01)
dd = 8:28
bsvol = sd(ret)
#
# reevaluate z with the code above
#
contour( dd , w0 , z , nlevels = 50 , zlim=c(9500,10500) )
```

```contour( dd , w0 , z , nlevels = 100 , zlim=c(9800,10000) )
```

```# -> d = 14  and w0 around 0.15:
max(z)                                        # 9873.06
logL( bsvol=sd(ret) , w0=0.15 , d=14 )        # 9872.65

# some more pictures:
image( dd , w0 , z )
```

```persp( dd , w0 , z )
```

```persp3D( dd , w0 , z )
```

```
# finally we check for bsvol:
w0 = seq(from=0.01,to=0.4,by=0.01)
bsvol = seq(from=0.1/100,to=4/100,by=0.1/100)
nbs = length(bsvol)
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] , d=14 )
}
}
contour( bsvol , w0 , z , nlevels = 50)
```

```contour( bsvol , w0 , z , nlevels = 80 , zlim=c(9800,10000))
```

```max(z)                                        # 9872.76
logL( bsvol=sd(ret) , w0=0.15 , d=14 )        # 9872.65

# -> Maximum at w0=15%, bsvol=sd(ret)=1.39%, d=14.

-----------------------------------------
#      ARCH(d)-Fit to SPX-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")   # returns look ok
```

```# compute spx-logL and make setup for contour-plot:
# first, we fix bsvol to its natural value:
bsvol = sd(ret)
bsvol
# we make a contour-plot in the (w0,d)-plane:
w0 = seq( from=0.05 , to=0.95 , by=0.05 )
dd = 1:40
nw = length(w0)
nd = length(dd)
z = matrix( 0 , nrow=nd , ncol=nw )
for(i in 1:nd)
{
for(j in 1:nw )
{
z[i,j] = logL( bsvol , w0[j] , dd[i] )
}
}
contour( dd , w0 , z , nlevels = 50 )
```

```contour( dd , w0 , z , nlevels = 80 , zlim=c(70000,72000) )
```

```# apparently w0 less than 0.3 and d between 8 and 28:
w0 = seq( from=0.01 , to= 0.3 , by=0.005 )
dd = 8:28
#
# reevaluate z = logL with the code above
#
contour( dd , w0 , z , nlevels = 100 , zlim=c(70000,72000) )
```

```contour( dd , w0 , z , nlevels = 200 , zlim=c(71000,72000) )
```

```# -> d = 14  and w0 around 0.15
max(z)                                         # 71413.99
logL( bsvol=sd(ret) , w0=0.15 , d=14 )         # 71413.54

# some more pictures:
persp( dd , w0 , z )
```

```persp3D( dd , w0 , z )
```

```# finally we check for bsvol:
w0 = seq( from=0.01 , to=0.3 , by=0.005 )
bsvol = seq( from=0.1/100 , to=2/100 , by=0.05/100 )
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] , d=14 )
}
}
contour( bsvol , w0 , z , nlevels = 100 , zlim=c(70000,72000) )
```

```contour( bsvol , w0 , z , nlevels = 200 , zlim=c(71000,72000) )
```

```# thus: Maximum at approximately
logL( bsvol=sd(ret) , w0=0.15 , d=14 )         # 71413.54
max(z)                                         # okay..
sd(ret)

-----------------------------------------
#       ARCH(d)-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")   # returns look ok
```

```# compute GE-logL and make setup for contour-plot:
# first, we fix bsvol to its natural value:
bsvol = sd(ret)
bsvol
# we make a contour-plot in the (w0,d)-plane:
w0 = seq( from=0.05 , to=0.95 , by=0.05 )
dd = 1:40
nw = length(w0)
nd = length(dd)
z = matrix( 0 , nrow=nd , ncol=nw )
for(i in 1:nd)
{
for(j in 1:nw )
{
z[i,j] = logL( bsvol , w0[j] , dd[i] )
}
}
contour( dd , w0 , z , nlevels = 50 )
```

```contour( dd , w0 , z , nlevels = 70 , zlim=c(50000,51000) )
```

```# apparently w0 less than 0.3 and d between 10 and 50:
w0 = seq( from=0.01 , to= 0.3 , by=0.005 )
dd = 10:50
#
# reevaluate z = logL with the code above
#
contour( dd , w0 , z , nlevels = 70 , zlim=c(50000,51000) )
```

```contour( dd , w0 , z , nlevels = 100 , zlim=c(50500,51000) )
```

```# -> d = 29  and w0 around 0.15
max(z)                                        # 50936.45
logL( bsvol=sd(ret) , w0=0.15 , d=29 )        # 50936.28

# some more pictures:
persp( dd , w0 , z )
```

```persp3D( dd , w0 , z )
```

```# finally we check for bsvol:
w0 = seq( from=0.01 , to=0.3 , by=0.005 )
bsvol = seq( from=0.5/100 , to=2.5/100 , by=0.05/100 )
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] , d=29 )
}
}
contour( bsvol , w0 , z , nlevels = 50)
```

```contour( bsvol , w0 , z , nlevels = 100 , zlim=c(50500,51000) )
```

```# thus: Maximum at approximately
logL( bsvol=sd(ret) , w0=0.15 , d=29 )        # 50936.28
max(z)                                        # 50936.61
sd(ret)   # 1.62%

```