It’s about fitting affine processes for interest rates. We are using one CIR-type process and two Vasiceks, which are correlated. These 3 processes evolve slowly over time stochastically. The evolution distribution provides the prior for the t+1 fitted value from the t one. The fitted interest rate at each time for each maturity is a linear function of the 3 processes, where the coefficients C(tau) and D_j(tau) vary by tau but are constant over time. They are fairly complex functions of the risk-neutal (rn) parameters for the processes. The code starts by computing those C and Ds. Then the real-world parameters are computed by adding in risk transforms lambda and psi to the rn parameters, and the processes are estimated at each time point t. Then the linear functions are computed.

We are not Stan experts and the code is pretty primitive. I’m sure it could be improved. There are known rotations of the parameters that can produce the same fits. We try to control that by using a less general model and choice of priors - empirically avoiding what appear to be alternate modes. The two Vasicek processes can be switched and our only control for that is to force one of them to have a higher variance. Anyway, code for our most general model with closed-form C and D functions is below. More general models need solutions of ODEs for C and D. We have code for that too.

data {

int<lower=0> N; //# rows

int<lower=0> U; //# columns

matrix[N,U] y;

real ts[U];

real opy; // observations per year

}

transformed data {

int z[U*N,2]; //index for rows and columns of each observation

int o; //# observations in y

o = 0;

for (n in 1:N) {

for (t in 1:U) {

o = o+1;

z[o,1] = n;

z[o,2] = t;

}

}

}

parameters {

real<lower=0.005, upper=0.8> kaprn1;

real<lower=0.005, upper=0.08> kaprn2;

real<lower=0.02, upper=0.0525> kaprn3;

real<lower=0.5, upper=2.2> omrn1;

real<lower=-40, upper=-7> omrn2;

real<lower=5, upper=45> omrn3;

real<lower=-0.6, upper=0.6> lb;

real<lower=-0.4, upper=0.75> ls2;

real<lower=0.02, upper=0.33> ls3; //diff sig3 - sig2

real<lower=-0.98, upper=-0.65> corr;

real<lower=-6, upper=-1> log_sigma_y; //log stdev of residuals

real<lower=-1.75, upper=-0.7> lam1;

real<lower=11, upper=40> lam2;

real<lower=-60, upper=-20> lam3;

real<lower=-1.25, upper=0.3> psi2;

real<lower=-3.5, upper=0> psi3;

real<lower=-6, upper=-0.5> ldel; //log of delta_0

real gam1;

vector<lower=0, upper=20>[N] rC; //CIR short rate

matrix<lower=-25, upper=25>[N,2] rV; //Vasicek short rate for each observed time point}

}

transformed parameters {

real beta;

vector[3] kap;

vector[3] om;

vector[3] sig;

vector[3] kaprn; //rn for risk-neutral

vector[3] omrn;

vector[U] tau;

vector[3] lam;

matrix[3,U] C; //C(T)

matrix[3,U] D; //D(T)

vector[U] W;

real del;

real sigma_y;

real h;

vector[U] corradj;

vector[2] Cmv; //first CIR mean and variance

vector[N-1] Cmean;

vector[N-1] Cvar;

matrix[N-1,2] Vmean; // means after first one for Vasicek processes

vector[2] Vm; //first Vasicek mean

vector[2] Vvar; //first Vasicek variance

matrix[2,2] Sig; // covariance matrix for Vasicek

matrix[N,U] meen; // fitted value at each point

//Compute C and D functions

del = exp(ldel);

beta = exp(lb);

kaprn[1] = kaprn1;

omrn[1] = omrn1;

sig[1] = 1;

kaprn[2] = kaprn2;

omrn[2] = omrn2;

sig[2] = exp(ls2);

kaprn[3] = kaprn3;

omrn[3] = omrn3;

sig[3] = sig[2]+ls3;

for (t in 1:U) tau[t] = ts[t];

h = (kaprn[1]^2 + 2*beta)^0.5;*

//W = 1/Q for CIR

for (t in 1:U) {

W[t] = 2h-(kaprn1+h)*(1-exp(h*tau[t])); //tau[t] is tau

D[1,t] = -2*(1-exp(h*tau[t]))/W[t]/tau[t]+gam1; //CIR*

D[2,t] = (1-exp(-kaprn2tau[t]))/kaprn2/tau[t]; //Vasicek

D[3,t] = (1-exp(-kaprn3*tau[t]))/kaprn3/tau[t]; //Vasicek*

}

for (t in 1:U) {

C[1,t] = -(omrn1/beta/tau[t])(2*log(2*h/W[t])+tau[t]*kaprn1+tau[t]**h) + del;*

C[2,t] = sig[2]^2(-2+2D[2,t]+kaprn[2]*tau[t]**D[2,t]^2)/4/kaprn[2]^2+omrn[2]/kaprn[2]*(1-D[2,t]);

C[3,t] = sig[3]^2(-2+2*D[3,t]+kaprn[3]*tau[t]**D[3,t]^2)/4/kaprn[3]^2+omrn[3]/kaprn[3]*(1-D[3,t]);

}

for (t in 1:U) corradj[t] = (D[2,t]+D[3,t]-1+(exp(-tau[t](kaprn[2]+kaprn[3]))-1)/(kaprn[2]+kaprn[3])/tau[t])*corr*sig[2]*sig[3]/kaprn[2]/kaprn[3];

// Compute short rate processes

lam[1] = lam1;

kap[1] = kaprn1 - beta * lam1;

om[1] = omrn1;

lam[2] = lam2;

om[2] = omrn2 + sig[2] * lam2;

kap[2] = kaprn[2]-sig[2]*psi2 ;

lam[3] = lam3;

om[3] = omrn3 + sig[3] * lam3;

kap[3] = kaprn[3]-sig[3]*psi3;

Cmv[1] = om[1]/kap[1];

Cmv[2] = om[1]/kap[1]^2*beta/2; //variance om1/kap1^2*beta/2

for (n in 1:N-1) {

Cmean[n] = rC[n]/exp(kap[1]/opy) + om[1]*(1-exp(-kap[1]/opy))/kap[1]; // mean for CIR process*

Cvar[n] = beta/kap[1](1-exp(-kap[1]/opy))*(rC[n]/exp(kap[1]/opy) + om[1]*(1-exp(-kap[1]/opy))/2)/kap[1];

} //gamma mean and variance for next CIR step, used in prior for next rC where step d = 1/opy

Vm[1] = om[2]/kap[2];

Vm[2] = om[3]/kap[3];

Vvar[1] = sig[2]^2/2/kap[2];

Vvar[2] = sig[3]^2/2/kap[3];

for (n in 1:N-1){ Vmean[n,1] = rV[n,1]+(om[2]-kap[2]*rV[n,1])/opy;

Vmean[n,2] = rV[n,2]+(om[3]-kap[3]*rV[n,2])/opy;} //mean for rV[n+1,i]*

for (i in 1:2) Sig[i,i] = sig[i+1]^2/opy;

Sig[1,2] = corrsig[3]*sig[2]/opy;

Sig[2,1] = Sig[1,2];

// Finally moments of the fit

for (n in 1:N) {

for (u in 1:U) {

meen[n,u] = D[1,u]*rC[n]+C[1,u]+D[2,u]*rV[n,1]+C[2,u]+D[3,u]*rV[n,2]+C[3,u]+corradj[u];

}

}

sigma_y = exp(log_sigma_y); //makes prior unbiased

//Priors not specified are uniform over their defined ranges, like non-negative reals for sigma_y; transformed parameters are simulated but don’t have priors.

}

model {

gam1 ~ double_exponential(0,0.01);

rC[1] ~ gamma(2*om[1]/beta,2*kap[1]/beta); // long term distribution of short rate

for (n in 2:N) rC[n] ~ gamma(Cmean[n-1]^2/Cvar[n-1],Cmean[n-1]/Cvar[n-1]); //in Stan mean of gamma is alpha/beta, etc

rV[1,1] ~ normal(Vm[1],Vvar[1]^0.5);

rV[1,2] ~ normal(Vm[2],Vvar[2]^0.5);

for (n in 2:N) {rV[n,] ~ multi_normal(Vmean[n-1,],Sig);}

for (u in 1:U) {for (n in 1:N) y[n,u] ~ normal(meen[n,u],sigma_y);}

}

generated quantities {

vector[o] log_lik;

for (j in 1:o) log_lik[j] = normal_lpdf(y[z[j,1], z[j,2]] | meen[z[j,1], z[j,2]], sigma_y);

}