---------------------------------------------------
#                                                 #
#     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 = read.table("C:/Users/Admin/Courses/FinancialStatistics/DAX.txt",header=TRUE,sep=";")
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:
spx = read.table("C:/Users/Admin/Courses/FinancialStatistics/spx.txt",header=TRUE,sep=";")
head(spx)
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:
ge = read.table("C:/Users/Admin/Courses/FinancialStatistics/GE.txt",header=TRUE,sep=";")
head(ge)
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%



Option Pricing
Time Series Models















Hochschule RheinMain,
Applied Mathematics:


<<previous chapter       back to chapter overview       next chapter>>