[R-SIG-Finance] Cox, Ingersoll, Ross/Vasicek parameter
Pedro Silva
psilva at porto.ucp.pt
Fri Apr 27 15:47:28 CEST 2007
Hi,
I did have the same problem and was not able to use the SSPIR or the
sde packages for the CIR model (for Vasicek's the SSPIR worked fine for me)
because the sde package assumes that the volatility is constant and although the SSPIR allows the user to define a function to specify
a time-varying volatility, as far as I can tell there is no way to
make the volatility dependent on the value of the state variable, as assumed by the CIR model.
So, I've programmed myself a small function, (in the bottom
of this message) which I called kalmanf, to perform the kalman
recursions with a state-dependent volatility. See if that is
useful to you (any feed-back is also very welcomed).
Note however that kalmanf does not estimate any parameter by itself,
only takes a set of parameters as given, together with the output of
a state-space model, and returns the output likelihood and, if required,
some additional information. In order to estimate the model parameters
you will need to combine the likelihood computations given by kalmanf with some optimization procedure using R optimization tools such as the optim or nlminb R routines. I've been trying that approach with mixed success, mainly because the likelihood function is very irregular and with many
local maxima and to get some reasonably credible values for the maximum-likelihood parameters, the optimization needs to be started many times from
different starting points.
I'll be very interested in hearing the experiences of anybody else who tried to solve the same problem.
Best regards,
A. Pedro Duarte Silva,
Associate Professor,
Fac. Economia e Gestão.
Univ. Católica Portuguesa at Porto
PORTO PORTUGAL
Kalmanf <- function(F,H,G=NULL,Qt,Qtpar,R,Yt,ut=NULL,z0=NULL,P0=NULL,result="all") {
# Runs the Kalman filter for the model
# Yt = H Zt + R wt
# Zt = F Zt-1 + G ut + Qt et
# where Zt represents the unobserved state vector,
# and wt, et are independent gaussian white noise vectors
#
# Written by A. Pedro Duarte Silva
# Faculdade de Econmia e Gestao - Univ. Catolica Portuguesa / Porto
# April 2007
# Parameters:
#
# F -- transition matrix (constant)
# H -- measurement coefficients matrix (constant)
# G -- coefficients of exogenous interventions (constant)
# Qt -- function that returns variable Choleski decompositions of the
# var-cov matrix for the noise in the transition equation
# Qt -- Fixed parameters (other than Zt) passed to the Qt function
# R -- matrix with the Choleski decomposition of var-cov for the noise
# in the measurement equation
# Yt -- matrix with the system output
# ut -- matrix (or vector) with system inputs (interventions)
# z0 -- vector with initial values for the state vector
# set to 0 by default
# P0 -- matrix with initial values for var-cov of the state vector
# set to "a large (diag(100)) value", representing difuse priors, by default
# result -- nature and amount of output returned
# Value:
# if result="all" a list with four components
# like -- -1 times the log-likelihood
# liket -- a vector with the observation contributions to the log-likelihood
# pred -- an array with the conditional means of Yt at time t-1
# Ytcov -- an array with the conditional var-cov of Yt at time t-1
# if result="like" -1 times the log-likelihood
# Parameter initializations
if (is.matrix(Yt)) {
T <- nrow(Yt) # number of periods
p <- ncol(Yt) # dimension of the output
}
else {
T <- length(Yt)
p <- 1
}
if (is.matrix(F)) n <- nrow(F) # dimension of the state vector
else n <- 1
if (!is.matrix(H)) dim(H) <- c(p,n)
# Allocate memory
FYt <- matrix(nrow=T,ncol=p) # Yt one-step forecasts (mean)
SigY <- array(dim=c(T,p,p)) # Yt one-step forecasts (var-cov)
Sigw <- matrix(nrow=p,ncol=p) # Var-Cov of the residuals in the measurement equation
Rwt <- array(dim=p) # Yt residuals
At <- matrix(nrow=n,ncol=p) # Auxiliary matrix (= FPt H' SigYt^(-1) )
c <- matrix(nrow=T,ncol=n) # constants in the transition equations
Zt <- array(dim=n) # state vectors
FZt <- array(dim=n) # Zt one-step forecasts (mean)
Qtet <- array(dim=n) # Zt residuals
Pt <- matrix(nrow=n,ncol=n) # Var-cov of the state vectors
FPt <- matrix(nrow=n,ncol=n) # Zt one-step forecasts (var-cov)
Sige <- matrix(nrow=n,ncol=n) # Var-Cov of the residuals in the transition equation (at time t)
lt <- array(dim=T) # time t contribution to the log-lokelihood
# Further Initializations
if (is.null(G)) c[] <- 0
else c <- ut %*% t(G)
Sigw <- R %*% t(R)
loglikcnst <- -0.5*( p*log(2*pi) )
if (is.null(z0)) Zt[] <- 0
else Zt <- z0
if (is.null(P0)) Pt <- diag(100,n,n)
else Pt <- P0
# Perform the Kalman recursions
for (t in 1:T) {
Siget <- Qt(Qtpar,Zt) %*% t(Qt(Qtpar,Zt))
FZt <- F %*% Zt + c[t,]
FPt <- F %*% Pt %*% t(F) + Siget
FYt[t,] <- H %*% FZt
Rwt <- Yt[t,] - FYt[t,]
SigY[t,,] <- H %*% FPt %*% t(H) + Sigw
At <- FPt %*% t(solve(SigY[t,,],H))
Qtet <- At %*% Rwt
Zt <- FZt + Qtet
Pt <- FPt - At %*% SigY[t,,] %*% t(At)
lt[t] <- loglikcnst -0.5*( log(det(SigY[t,,])) + Rwt %*% solve(SigY[t,,],Rwt) )
}
# Compute the log-likelihood and return the results
logLik <- sum(lt)
if (result == "like") return(-logLik)
if (result == "all") return(list(like=-logLik,liket=lt,pred=FYt,Ytcov=SigY))
}
Message: 3
Date: Fri, 27 Apr 2007 00:02:24 +0200
From: "Arne Krombach" <arne.krombach at web.de>
Subject: [R-SIG-Finance] Cox, Ingersoll, Ross/Vasicek parameter
estimation via Kalman-Filter (SSPIR)
To: <r-sig-finance at stat.math.ethz.ch>
Message-ID: <E1HhC2F-000072-00 at smtp06.web.de>
Content-Type: text/plain
Dear R-Users,
I am trying to estimate the parameters for a CIR 1-/2-/3-Factor model via
Kalman filtering. The idea was to do it via the package SSPIR, but I have
problems to transform the CIR state space form into the SSPIR syntax.
Actually, I am not quite sure if this is possible at all.
Did anybody already realise a CIR/Vasicek -parameter estimation via R?
To get an idea how a CIR-1-factor state space form looks like (Y= vector of
zerobond-yields with the maturities tau at month t):
Y(t,tau)= ( -1/tau * log A(tau) ) + (1/tau * B(tau) * x(t) + e
Where :
A(tau) = [(2 * h * e^( a + lambda + h) * tau/2) / (2 * h + (a + lambda + h)
* (e^(tau * h) - 1)] ^(2 * b/ sigma^2)
B(tau) = [(2 * (e^(tau * h) -1)) / 2 * h + (a + lambda + h) * (e^(tau * h)
- 1)]
h = sqrt[(a + lambda)^2 + 2*sigma^2]
which means that, besides the time varying factor x, the parameters a, b,
sigma, and lambda have to be estimated via Kalman filtering.
Is it possible to implement this in SSPIR? If you have any other suggests or
if you already implement comparable models in R, I would be very grateful
for your help.
Kind Regards,
Arne Krombach
Esta mensagem (incluindo quaisquer anexos) pode conter informação confidencial ou legalmente protegida para uso exclusivo do destinatário. Se não for o destinatário pretendido da mesma, não deverá fazer uso, copiar, distribuir ou revelar o seu conteúdo (incluindo quaisquer anexos) a terceiros, sem a devida autorização. Se recebeu esta mensagem por engano, por favor informe o emissor, por e-mail e elimine-a imediatamente. Obrigado.
This message may contain confidential information or privile...{{dropped}}
More information about the R-SIG-Finance
mailing list