There is an obscure and useful function that makes it easy to fit stochastic differential equations to data insofar as the model can be linearized without causing too much trouble. The function discretizes the continuous-time (i.e., differential equation) state matrices A, the drift or state transition matrix, B, the input or covariate coefficient matrix, and Q, diffusion or noise covariance matrix. That means that the function essentially takes the differential equation in matrix form and solves it for a given time step. The discretized matrices function like those of an autoregressive process. Some details of this approach and what this does can be found here but not exactly a complete implementation, namely with matrix B. So this is one of those code blocks I just have backed up in several project folders in various languages.
Thanks to Mike Hunter of Georgia Tech for originally creating this and providing it in OpenMx.
R:
discretize <- function(A, B, Q, deltaT){
deltaT <- c(deltaT)
I <- diag(ncol(A))
BLOCK <- matrix(0, nrow=2*nrow(A), ncol=2*ncol(A))
# First Block expm for A integral, and expm(A*deltaT)
BLOCK[1:(2*nrow(A)), 1:ncol(A)] <- 0
BLOCK[1:nrow(A), (nrow(A)+1):(2*nrow(A))] <- I
BLOCK[(nrow(A)+1):(2*nrow(A)), (nrow(A)+1):(2*nrow(A))] <- A
BLOCK <- OpenMx::expm(BLOCK*deltaT)
expA <- BLOCK[(nrow(A)+1):(2*nrow(A)), (nrow(A)+1):(2*nrow(A))]
intA <- BLOCK[1:nrow(A), (nrow(A)+1):(2*nrow(A))]
# Second Block expm for discretized Q
BLOCK[1:(nrow(A)), 1:ncol(A)] <- -t(A)
BLOCK[(nrow(A)+1):(2*nrow(A)), 1:ncol(A)] <- 0
BLOCK[1:nrow(A), (nrow(A)+1):(2*nrow(A))] <- Q
BLOCK[(nrow(A)+1):(2*nrow(A)), (nrow(A)+1):(2*nrow(A))] <- A
BLOCK <- OpenMx::expm(BLOCK*deltaT)
Bd <- intA %*% B
Qd <- t(expA) %*% BLOCK[1:nrow(A), (nrow(A)+1):(2*nrow(A))]
return(list(Ad=expA, Bd=Bd, Qd=Qd))
}
C++ (RcppArmadillo):
// [[Rcpp::depends(RcppArmadillo)]]
#include <RcppArmadillo.h>
using namespace Rcpp;
using namespace arma;
using namespace std;
// [[Rcpp::export]]
List discretize(mat A, mat B, mat Q, double dt) {
mat M, expA, intA, Ad, Bd, Qd, I;
List output;
I.eye(A.n_rows, A.n_cols);
M.zeros(2*A.n_rows, 2*A.n_cols);
// First Block expm for A integral, and expm(A*deltaT)
M.submat(0, A.n_rows, A.n_rows-1, 2*A.n_rows-1) = I;
M.submat(A.n_rows, A.n_rows, 2*A.n_rows-1, 2*A.n_rows-1) = A;
M=expmat(M*dt);
intA = M.submat(0, A.n_rows, A.n_rows-1, 2*A.n_rows-1);
expA = M.submat(A.n_rows, A.n_rows, 2*A.n_rows-1, 2*A.n_rows-1);
// Second Block expm for discretized Q
M.submat(0, 0, A.n_rows-1, A.n_rows-1) = -1*A.t();
M.submat(A.n_rows, 0, 2*A.n_rows-1, A.n_rows-1) = zeros(A.n_rows, A.n_cols);
M.submat(0, A.n_rows, A.n_rows-1, 2*A.n_rows-1)=Q;
M.submat(A.n_rows, A.n_rows, 2*A.n_rows-1, 2*A.n_rows-1) = A;
M=expmat(M*dt);
output("Ad")=expA;
output("Bd")=intA*B;
output("Qd")=expA.t() * M.submat(0, A.n_rows, A.n_rows-1, 2*A.n_rows-1);
return output;
}
Stan:
functions{
matrix[] discretize(matrix A, matrix B, matrix Q, real dt) {
int m = rows(A);
int n = m;
matrix[m,n] expA;
matrix[m,n] intA;
matrix[m,n] I = diag_matrix(rep_vector(1.0, m));
matrix[m,n] output[3];
// First Block expm for A integral, and expm(A*deltaT)
matrix[2*m,2*n] M = append_col(rep_matrix(0, 2*m, n), append_row(I, A));
M = matrix_exp(M*dt);
intA = block(M, 1, n+1, m, n);
expA = block(M, m+1, n+1, m, n);
// Second Block expm for discretized Q
M = append_col( append_row(-1*A', rep_matrix(0, m, n)), append_row(Q, A));
M = matrix_exp(M*dt);
output[1] = expA;
output[2] = intA*B;
output[3] = expA' * block(M, 1, n+1, m, n);
return output;
}
}