' shm
' Calibration of Schwartz Smith model parameters

' path:
   %path = @runpath
   cd %path

' WF for Schwartz-Smith model
   wfcreate(wf=ssk) W 01/05/1990 02/17/1995

' Import data
   import futureslog.xlsx

' Initial values

   c(1) = 2 ' K
   c(2) = 0.2 ' Sigmax
   c(3) = 0.2 ' Lambdax
   c(4) = 0.2 ' Sigmae
   c(5) =  0.02 ' Mu
   c(6) = 0.02 ' Mu star
   c(7) =  0.2 ' Pxe
   c(8)= 0.02 ' forward1 vol
   c(9)= 0.02 ' forward2 vol
   c(10)= 0.02 ' forward3 vol
   c(11)= 0.02 ' forward4 vol
   c(12)= 0.02 ' forward5 vol
  

' List of series
   %Listofseries= @wlookup("*m*","series")

' Number of series
   scalar nf
   nf = @wcount(%Listofseries)

' Number of observations
   scalar n
   n = @obs(m1) 

' Maturities
   rowvector(nf) matur
   matur.fill 1/12, 5/12, 9/12, 13/12, 17/12 
   
' Time step between observations
   scalar dt
   dt=5/252

' Initial state variables and cov. matrix
   vector(2) svec0 
   svec0.fill 0, 3.1307  
   sym(2) svar0 
   svar0.fill 0.1, 0.05, 0.2 



' Sspace 
   sspace estimationss

' Set initial values in sspace
   estimationss.append @mprior svec0 
   estimationss.append @vprior svar0 


' Transition Equation. Xt=c+Q(t-1) + Nt. Eq. 27 Paper, page18. 

' Xt = X(t-1) * exp(-k*dt ) + var(Xt)  
   estimationss.append @state ShTDev= ShTDev(-1) * exp(-C(1) * (dt)) + [ename = e1]
' Et = E(t-1) + µ*dt  + var(Et)   
   estimationss.append @state EquilPrice= EquilPrice(-1) + C(5) * (dt) + [ename = e2]

' Variance of Xt = (1 - exp(-2*k*dt ) * (1/(2*k)) * Sigmax^2
   estimationss.append @evar var(e1) = (1-exp(-2*C(1)*(dt))) * 0.5 * (C(2) * C(2) / C(1))
' Variance of Et = Sigmae^2
   estimationss.append @evar var(e2) = C(4) * C(4) * (dt)
' Covariance Xt,Et = (1 - exp(-k*dt ) * Pxe * Sigmae * Sigmax / k
   estimationss.append @evar cov(e1,e2) = (1-exp(-C(1)*(dt))) * (C(7) * C(2) * C(4) / C(1))


' Measurement Equation  Yt=A(Ti) *Zt*Xt + Et. Eq. 28 Paper, page 19.
   estimationss.append @signal m1 = exp(-C(1)*(matur(1))) * ShTDev+ EquilPrice+ C(6) * (matur(1)) - (1-exp(-C(1)*(matur(1)))) * (C(3) / C(1)) + (1/2) * (1-exp(-2*C(1) * (matur(1)))) * ((C(2) * C(2)) / (2*C(1))) + (1/2) * (C(4) * C(4) * matur(1)) + (1-exp(-C(1)*(matur(1)))) * C(7) * C(2) * C(4) / C(1) + [ename = e3]

   estimationss.append @signal m5 = exp(-C(1)*(matur(2))) * ShTDev+ EquilPrice+ C(6) * (matur(2)) - (1-exp(-C(1)*(matur(2)))) * (C(3) / C(1)) + (1/2) * (1-exp(-2*C(1) * (matur(2)))) * ((C(2) * C(2)) / (2*C(1))) + (1/2) * (C(4) * C(4) * matur(2)) + (1-exp(-C(1)*(matur(2)))) * C(7) * C(2) * C(4) / C(1) + [ename = e4]

   estimationss.append @signal m9 = exp(-C(1)*(matur(3))) * ShTDev+ EquilPrice+ C(6) * (matur(3)) - (1-exp(-C(1)*(matur(3)))) * (C(3) / C(1)) + (1/2) * (1-exp(-2*C(1) * (matur(3)))) * ((C(2) * C(2)) / (2*C(1))) + (1/2) * (C(4) * C(4) * matur(3)) + (1-exp(-C(1)*(matur(3)))) * C(7) * C(2) * C(4) / C(1) + [ename = e5]

   estimationss.append @signal m13 = exp(-C(1)*(matur(4))) * ShTDev+ EquilPrice+ C(6) * (matur(4)) - (1-exp(-C(1)*(matur(4)))) * (C(3) / C(1)) + (1/2) * (1-exp(-2*C(1) * (matur(4)))) * ((C(2) * C(2)) / (2*C(1))) + (1/2) * (C(4) * C(4) * matur(4)) + (1-exp(-C(1)*(matur(4)))) * C(7) * C(2) * C(4) / C(1) + [ename = e6]

   estimationss.append @signal m17 = exp(-C(1)*(matur(5))) * ShTDev+ EquilPrice+ C(6) * (matur(5)) - (1-exp(-C(1)*(matur(5)))) * (C(3) / C(1)) + (1/2) * (1-exp(-2*C(1) * (matur(5)))) * ((C(2) * C(2)) / (2*C(1))) + (1/2) * (C(4) * C(4) * matur(5)) + (1-exp(-C(1)*(matur(5)))) * C(7) * C(2) * C(4) / C(1) + [ename = e7]

' Forwards volatilities
   estimationss.append @evar var(e3) = C(8) ^ 2
   estimationss.append @evar var(e4) = C(9) ^ 2
   estimationss.append @evar var(e5) = C(10) ^ 2
   estimationss.append @evar var(e6) = C(11) ^ 2
   estimationss.append @evar var(e7) = C(12) ^ 2

   estimationss.ml


' Show results
   show estimationss