SsfPack S+FinMetrics - Solving parameters for the Vasicek Model

["Ellerton, Tim" <Tim.Ellerton@morganstanley.com>]

I am trying to use a Kalman filter/ MLE to solve the parameters for a single-factor Vasicek model and the sigma seems to be consistently out. 

 

The standard OLS and MLE methods are consistently more accurate.

 

Any help/guidance would be much appreciated.

 

Many thanks,

 

 

Tim

 

 

# Title: Kalman/MLE Estimation

#Load modules

module(finmetrics)

vasicek.ssf = function(param, tau=NULL, freq=1/252){

#Zivot, Wang and Koopman (2003) - "State Space Modeling in Macroeconomics and Finance Using SsfPack for S+Finmetrics"

## 1. Check for valid input.

if (length(param) < 5)

stop("param must have length greater than 4.")

N = length(param) - 4

if (length(tau) != N)

stop("Length of tau is inconsistent with param.")

## 2. Extract parameters and impose constraints.

Kappa = exp(param[1]) ## Kappa > 0

Theta = param[2]

Sigma = exp(param[3]) ## Sigma > 0

Lamda = param[4]

Var = exp(param[1:N+4]) ## meas eqn stdevs

## 3. Compute Gamma, A, and B.

Gamma = Theta + Sigma * Lamda / Kappa - Sigma^2 / (2 * Kappa^2)

B = (1 - exp(-Kappa * tau)) / Kappa

lnA = Gamma * (B - tau) - Sigma^2 * B^2 / (4 * Kappa)

## 4. Compute a, b, and Phi.

a = Theta * (1 - exp(-Kappa * freq))

b = exp(-Kappa * freq)

Phi = (Sigma^2 / (2 * Kappa)) * (1 - exp(-2 * Kappa * freq))

## 5. Compute the state space form.

mDelta = matrix(c(a, -lnA/tau), ncol=1)

mPhi = matrix(c(b, B/tau), ncol=1)

mOmega = diag(c(Phi, Var^2))

## 6. Duan and Simonato used this initial setting.

A0 = Theta

P0 = Sigma * Sigma / (2*Kappa)

mSigma = matrix(c(P0, A0), ncol=1)

## 7. Return state space form.

ssf.mod = list(mDelta=mDelta, mPhi=mPhi, mOmega=mOmega, mSigma=mSigma)

CheckSsf(ssf.mod)

}

fit.vasicek.kalman.mle <- function(x, est.kappa = 5, est.theta = 0, est.sigma = 0.1, est.lamda = 0, est.spread = 0.01, tau = 1, freq = 1/252){

#calculate MLE's of unknown vasicek parameters in state space model

init.vasicek <- c(log(est.kappa), est.theta, log(est.sigma), est.lamda, log(est.spread))

rowIds(init.vasicek) <- c("ln.kappa", "theta", "ln.sigma", "lamda", "ln.spread")

vasicek <- SsfFit(parm = init.vasicek, data = "" FUN=vasicek.ssf, tau=tau, freq=freq, trace=F, control=nlminb.control(abs.tol=1e-6, rel.tol=1e-6, x.tol=1e-6, eval.max=1000, iter.max=500))

#summary(vasicek)

ssf.fit = vasicek.ssf(vasicek$parameters, tau = est.tau, freq = freq)

auto.cor.test <- autocorTest(KalmanFil(x, ssf.fit)$std.innov, method = "lb")

#plot(x)

#plot(SsfCondDens(x, ssf.fit, task="STSMO", ikf=NULL))

#plot(KalmanFil(x, ssf.fit))

result <- c(exp(vasicek$parameters["ln.kappa"]), vasicek$parameters["theta"], exp(vasicek$parameters["ln.sigma"]), vasicek$parameters["lamda"], exp(vasicek$parameters["ln.spread"]), auto.cor.test$p.value)

rowIds(result) <- c("kappa", "theta", "sigma", "lamda", "spread", "ljung.box.p")

return(result)

}

fit.vasicek.ols <- function(x, freq = 1/252){

#calculate unknown vasicek parameters by linear regression

#http://www.sitmo.com/doc/Calibrating_the_Ornstein-Uhlenbeck_model

x <- seriesMerge(x, tslag(x, trim=T))

colIds(x) <- c("S","S.lag")

x.ols <- OLS(S~S.lag, data = "">

a <- x.ols$coef["S.lag"]

b <- x.ols$coef["(Intercept)"]

e <- summary(x.ols)$sigma

kappa <- - log(a) / freq

theta <- b / (1 - a)

sigma <- e * sqrt((-2 * log(a)) / (freq * (1 - a * a)))

result <- c(kappa, theta, sigma)

rowIds(result) <- c("kappa", "theta", "sigma")

return(result)

}

 

fit.vasicek.mle <- function(x, day.count = 252){

#calculate unknown vasicek parameters using mle

#http://www.sitmo.com/doc/Calibrating_the_Ornstein-Uhlenbeck_model

x <- seriesMerge(x, tslag(x, trim=T))

colIds(x) <- c("S","S.lag")

n <- nrow(x)

dt <- 1 / day.count

Sx <- sum(x[, "S.lag"])

Sxx <- sum(x[, "S.lag"]^2)

Sy <- sum(x[, "S"])

Sxy <- sum(x[, "S.lag"] * x[, "S"])

Syy <- sum(x[, "S"]^2)

theta <- (Sy*Sxx - Sx*Sxy) / ( n*(Sxx - Sxy) - (Sx^2 - Sx*Sy))

kappa <- -log( (Sxy - theta * Sx - theta * Sy + n * theta^2) / (Sxx - 2 * theta * Sx + n * theta^2) ) / dt

a <- exp(-kappa * dt)

sigma.h.2 <- (Syy - 2 * a * Sxy + a^2 * Sxx - 2 * theta * (1 - a) * (Sy - a * Sx) + n * theta^2 * (1 - a)^2) / n

sigma = sqrt(sigma.h.2 * 2 * kappa / (1 - a^2))

result <- c(kappa, theta, sigma)

rowIds(result) <- c("kappa", "theta", "sigma")

return(result)

}

#simulate data

Positions <- timeSeq(from=timeDate("01/01/2008"), to=timeDate("01/01/2009"), by="days", k.by=1)

kappa <- 15

theta <- 0

sigma <- 0.12

x <- OU.gensim(rho = c(kappa, theta, sigma), n.sim = length(Positions), n.burn = 500, aux = OU.aux(X0 = 0.1, ndt = 1, t.per.sim = 1/252))

x <- timeSeries(x,Positions)

#plot(x)

#estimate parameters

vasicek.ols <- fit.vasicek.ols(x)

vasicek.mle <- fit.vasicek.mle(x)

vasicek.kalman.mle <- fit.vasicek.kalman.mle(x, vasicek.mle["kappa"], vasicek.mle["theta"], vasicek.mle["sigma"], freq = 1/252)

#print parameters

vasicek.ols

vasicek.mle

vasicek.kalman.mle

---

NOTICE: If received in error, please destroy, and notify sender. Sender does not intend to waive confidentiality or privilege. Use of this email is prohibited when received in error. We may monitor and store emails to the extent permitted by applicable law.