Here is one of the codes.....
Thanks a lot...
module(finmetrics)
genvasicek.ssf = function(param, tau=NULL, freq=1/12){
## 1. Check for valid inputs
if(length(param) < 5)
stop("Parameters must have length greater than 4")
N = length(param) - 13
if (length(tau) != N)
stop("Length of Tau is inconsistent with Parameters")
## 2.Extract Parameters and Impose Constraints
k1 = (param[1])
k2 = (param[2])
k3 = (param[3])
lambda1 = param[4]
lambda2 = param[5]
lambda3 = param[6]
Delta = (param[7])
s1 = (param[8])
s2 = (param[9])
s3 = (param[10])
r21 = param[11]
r31 = param[12]
r32 = param[13]
st1 = (param[14])
st2 = (param[15])
st3 = (param[16])
st4 = (param[17])
## 3. Compute Matrixes
A1 = function(x,tau){
(1-exp(-x*tau))/x
}
u1 = -A1(k1,tau) ## ok
u2 = -A1(k2,tau) ## ok
u3 = -A1(k3,tau) ## ok
H = cbind (-u1/tau,-u2/tau,-u3/tau) ## ok
v1 = (lambda1/k1)*(tau - A1(k1,tau)) + (lambda2/k2)*(tau - A1(k2,tau)) +
(lambda3/k3)*(tau - A1(k3,tau))
v2 = +3 * Delta * tau # Teste
v311 = ((s1*s1)/(k1*k1)) * (tau - A1(k1,tau) - A1(k1,tau) +
A1(2*k1,tau))
v312 = ((s1*s2*r21)/(k1*k2)) * (tau - A1(k1,tau) - A1(k2,tau) +
A1(k1+k2,tau))
v313 = ((s1*s3*r31)/(k1*k3)) * (tau - A1(k1,tau) - A1(k3,tau) +
A1(k1+k3,tau))
v322 = ((s2*s2)/(k2*k2)) * (tau - A1(k2,tau) - A1(k2,tau) +
A1(2*k2,tau))
v321 = ((s1*s2*r21)/(k1*k2)) * (tau - A1(k1,tau) - A1(k2,tau) +
A1(k1+k2,tau))
v323 = ((s2*s3*r32)/(k2*k3)) * (tau - A1(k2,tau) - A1(k3,tau) +
A1(k2+k3,tau))
v333 = ((s3*s3)/(k3*k3)) * (tau - A1(k3,tau) - A1(k3,tau) +
A1(2*k3,tau))
v331 = ((s1*s3*r31)/(k1*k3)) * (tau - A1(k3,tau) - A1(k1,tau) +
A1(k3+k1,tau))
v323 = ((s2*s3*r32)/(k2*k3)) * (tau - A1(k2,tau) - A1(k3,tau) +
A1(k2+k3,tau))
v = v1+v2+0.5*(v311+v312+v313+v322+v321+v323+v333+v331+v323)
d = -v
A = matrix(c(0,0,0,0,0,0,0,0,0),ncol=3)
A[1,1] = 1-k1*freq
A[1,2] = 0
A[1,3] = 0
A[2,1] = 0
A[2,2] = 1-k2*freq
A[2,3] = 0
A[3,1] = 0
A[3,2] = 0
A[3,3] = 1-k3*freq
ct = c(-lambda1*freq,-lambda2*freq,-lambda3*freq)
Q = matrix(c(0,0,0,0,0,0,0,0,0),ncol=3)
Q[1,1] = s1^2
Q[1,2] = s1*s2*r21
Q[1,3] = s1*s3*r31
Q[2,1] = s2*s1*r21
Q[2,2] = s2^2
Q[2,3] = s2*s3*r32
Q[3,1] = s1*s3*r31
Q[3,2] = s2*s3*r32
Q[3,3] = s3^2
## 4. Compute the State Space Form
mDelta = matrix(c(ct,d),ncol=1)
mPhi = rbind(A,H)
mOmega1 = matrix(c(Q,0,0,0,0,0,0,0,0,0,0,0,0),ncol=7)
mOmega2 matrix(c(0,0,0,0,0,0,0,0,0,0,0,0,diag(c(st1,st2,st3,st4))),ncol=7)
mOmega = rbind(mOmega1,mOmega2)
mSigma = matrix(c(0,0,0,0,0,0,0,0,0,0,0,0),ncol=3)
mSigma[1,1] = ((s1^2)/(2*k1))*(1-exp(-2*k1*freq))
mSigma[2,2] = ((s2^2)/(2*k2))*(1-exp(-2*k2*freq))
mSigma[3,3] = ((s3^2)/(2*k3))*(1-exp(-2*k3*freq))
## 7. Return State Space Form
ssf.mod = list(mDelta=mDelta, mPhi=mPhi, mOmega=mOmega)
CheckSsf(ssf.mod)
}
start = c(0.1,0.1,0.1,0.3,0.3,0.3,0.06,0.005,0.005,0.005,0.001,0.001,0.001,
0.003,0.1,0.003,0.01)
names(start)=c("logkappa1","logkappa2","logkappa3","lambda1","lambda2","lambda3","logDelta","logSigma1","logSigma2","logSigma3","Corr12","Corr13","Corr23","logSigmat1","logSigmat2","logSigmat3","logSigmat4")
tau = c(0.25,0.5,1,5)
l = c(0,0,0,0,0,0,0,0,0,-0,-1,-1,-1,0,0,0,0)
u = c(10,10,10,10,10,10,0.5,10,10,10,1,1,1,10,10,10,10)
ans.vasicek = SsfFit(start, fama.bliss, genvasicek.ssf, tau=tau, freq=1/12,
trace=T, lower=l, upper=u, control=nlminb.control(abs.tol=1e-6,,
rel.tol=1e-6, x.tol=1e-6, eval.max=1000,iter.max=135))
On 8/13/07, Bernardo Ribeiro <ribeirob@gmail.com>
wrote:>
> Hey all,
>
> I am trying to work under a State Space form, but I didn't get the help
> exactly.
> Have anyone eles used this functions?
>
> I was used to work with S-PLUS, but I have some codes I need to adpt.
>
> Thanks alot,
>
> Bernardo
>
[[alternative HTML version deleted]]