R Labs 7

Optimization heuristics in finance

R Example 7.1 (Estimating GARCH(1,1) with Simulated Annealing): In R the function  optim includes an option for performing optimization based on simulated annealing. To use it, set method="SANN". Then define the function to be minimized (or maximized), and pass this function to  optim,   together with a vector of initial values for parameters to be optimized over,  and appropriate control settings (i.e. the type of optimization -if a maximization or a minimization-, the maximum number of iterations, etc.).

Recall that the GARCH(1,1) model for a log-return time series \{r_t\}_{t=1}^T has the form

\begin{eqnarray*}
r_t &=& \mu + a_t,\ \quad a_t \sim N(0,\sigma_t^2) \\
\sigma^2_t &=& \omega + \alpha_1\cdot a_{t-1}^2 + \beta_1 \cdot \sigma_{t-1}^2
\end{eqnarray*}
with -1\le \mu \le 1,\, \omega > 0,\, \alpha_1 \ge 0,\, \beta_1 \ge 0,\, \alpha_1 + \beta_1 < 1. Given the above constraints we want the values for \psi = (\omega, \alpha_1, \beta_1, \mu) that maximizes
the log likelihood function (without additive constants)
\begin{equation*}
F(\psi) = -\frac{1}{2} \sum_{t=1}^T \left(\ln \sigma_t^2 + \frac{a_t^2}{\sigma_t^2}\right)
\end{equation*}

The R code for maximizing F(\psi) with Simulated Annealing follows. The results are compared to results obtained by the  garchFit method.

###  @Argimiro Arratia, 2013
### Maximum Likelihood methods for GARCH(1,1), and SANN (Simulated Annealing)
######
##psi = c(a0,a1,b1,mu) #coef to maximize
## r = c(r1,r2,..., rn) ##log return series
## the objective function (This is the solution to Exercise 7.7.3 of the book)
fn   a0 = psi[1]
  a1 = psi[2]
  b1 = psi[3]
  mu = psi[4] ##mean(r) ##In case we want to omit mu as parameter
  n = length(r)
  v = rep(0,n)
  A = rep(0,n)
  v[1] = r[1]^2
  A[1] = r[1]-mu
  sum = log(v[1]) + A[1]^2/v[1]
  for (i in 2:n){
    A[i] = r[i]-mu
    v[i] = a0 + a1*A[i-1]^2 + b1*v[i-1]
    sum = sum + log(v[i]) + (A[i]^2/v[i])
    }
  return(-sum)
}

##Set up an indicator function for the GARCH(1,1) constraints for parameters
I  0 && psi[2]>= 0 && psi[3] >= 0 && psi[2]+psi[3] < 1      && psi[4]>-1/2 && psi[4]<1/2) ##psi[4] includes mu as parameter
  {return(1)}
  else {return(0)}
}

LI   return(fn(psi,r)*I(psi))
}
##follow optimization method, use with specific r and initial parameters parinit
#set working directory where you have saved file ibexRlog (see Data)
setwd("your-working-directory")

##load the Ibex file with log returns for SAN, TEF, REP, BBVA
ibexRl = read.table("ibexRlog", header=T)

sanRl <-ibexRl$sanR ##extract log returns of SAN
##remove undefined terms NA
sanRl r = rev(sanRl) ##reverse sanRl since begins with last day and to the past
mu = mean(r)
parinit = c(0.1,0.1,0.1,0.1) ##,mu)
opt = optim(parinit,LI,r=r,method = "SANN", control= list(fnscale=-1,maxit=5000))
opt$par
opt$value ##returns the value of fn corresponding to par

##To know the current (or day n) estimated variance, we must know variance on day n-1:
# n = length(r)
# sigma2[n]= opt$par[1] + opt$par[2]*(r[n-1]-mu)^2 + opt$par[3]*sigma2[n-1]
##To annualize this estimated today's volatility
# sigan= sqrt(252)*sqrt(sigma2[n])
##Instead we are contempt with a mean (or long-run average) volatility per day:
Vl = sqrt(opt$par[1]/(1-opt$par[2]-opt$par[3]))

##Results for r = rev(sanRl) with
##i) without mu as parameter i.e mu = 0 (comment out mu in fn and cut parinit to 3 first parameters, Uncomment 3 lines below)
# opt$par
# Vl = sqrt(opt$par[1]/(1-opt$par[2]-opt$par[3]))
# Vl

## ii) mu as parameter (uncomment 4 lines below)
# opt$par
# mu
# Vl = sqrt(opt$par[1]/(1-opt$par[2]-opt$par[3]))
# Vl

##iii) with mu = mean(r) fixed (not a parameter so is not updated ever)
opt$par
Vl = sqrt(opt$par[1]/(1-opt$par[2]-opt$par[3]))
Vl

##Warning: have in mind that another run will give different coef. due to random nature of SANN method
##The important thing is that Vl has value in (.025, 0.03)
##These seem to be better attained with mu= mean(r) (constant) or mu as parameter
##All experiments above are compared to  built in garch function. Remember to load
library(fBasics)
library(fGarch)
m2=garchFit(~garch(1,1),data=r,trace=F,cond.dist="std")
summary(m2)
##Results om, alfa1, beta1. Compute volatility Vll
Vll = sqrt(om/(1-alfa1-beta1))