[R-SIG-Finance] Cox, Ingersoll, Ross/Vasicek parameter

[Pedro Silva]

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