The following page documents the simulation happening in the model described here, i.e. the Wishart based stochastic volatility. The parameters for the Euler-Maruyama simulation are:

Example of Arbitrage free volatility surface (dynamics of vol surface) generated by the model (here x axis is tenor structure, y axis is strike level and z axis call option price):

The simulation of Wishart varinace process is accomplished via the following code:

function [S,eigS,timestep,minh]=Wishart(T,h,d,alpha,Q,M,s_0)
%
% Simulates the d-dimensional Wishart process on the interval [0,T]
% that follows the stochastic differential equation
% dS_t=sqrt{S_t}*dB_t*Q+Q'*dB_t'*sqrt{S_t}+(S_t*M+M'*S_t+alpha*Q'*Q)dt
% with initial condition S_0=s_0.
%
% Method of discretization: Euler-Maruyama
% Starting step size: h
% In order to guarantee positive semidefiniteness of S, the step size will be
% reduced iteratively if necessary.
%
% Output:
% S is a three dimensional array of the discretized Wishart process.
% eigS is a matrix consisting the eigenvalues of S
% timestep is the vector of all timesteps in [0,T]
% minh is the smallest step size used, i.e. minh=min(diff(timestep))
%--------------------------------------------------------------------
horg=h;
minh=h;
timestep=0;
[V_0,D_0]=eig(s_0);
eigS=sort(diag(D_0)');
% to ensure evolution in the right space
drift_fix=alpha*Q'*Q;
B_old=0;
B_inc=sqrt(h)*normrnd(0,1,d,d);
% not clear why he needs Q : the correlation structure?
vola=V_0*sqrt(D_0)*V_0'*B_inc*Q;
drift=s_0*M;
S_new = s_0+vola+vola'+(drift+drift'+drift_fix)*h;
[V_new,D_new]=eig(S_new);
eigS=[eigS;sort(diag(D_new)')];
% concatenate along third dimension
S=cat(3,s_0,S_new);
t=h;
timestep=[timestep;t];
flag=0;1
while t+h<T,
B_old=B_old+B_inc;
B_inc=sqrt(h)*normrnd(0,1,d,d);
S_t=S_new; V_t=V_new; D_t=D_new;
sqrtm_S_t=V_t*sqrt(D_t)*V_t';
vola=sqrtm_S_t*B_inc*Q;
%vola= V_t*sqrt(D_t)*V_t'*B_inc*Q;
drift=S_t*M;
S_new = S_t+vola+vola'+(drift+drift'+drift_fix)*h;
[V_new,D_new]=eig(S_new);
mineig=min(diag(D_new));
% check whether the minimum eigenvalue is negative
while mineig<0,
h=h/2;
minh=min(minh,h);
flag=1;
if h<eps, error('Step size converges to zero'), return, end
B_inc=0.5*B_inc+sqrt(h/2)*normrnd(0,1,d,d);
vola=sqrtm_S_t*B_inc*Q;
drift=S_t*M;
S_new = S_t+vola+vola'+(drift+drift'+drift_fix)*h;
mineig=min(eig(S_new));
end
if flag==0,
eigS=[eigS;sort(diag(D_new)')];
S=cat(3,S,S_new);
t=t+h;
timestep=[timestep;t];
end
if flag==1,
[V_new,D_new]=eig(S_new);
eigS=[eigS;sort(diag(D_new)')];
S=cat(3,S,S_new);
t=t+h;
timestep=[timestep;t];
flag=0;
h=horg;
end
end
h_end=T-t;
if h_end>0,
B_inc=sqrt(h_end)*normrnd(0,1,d,d);
S_t=S_new; V_t=V_new; D_t=D_new;
vola=V_t*sqrt(D_t)*V_t'*B_inc*Q;
drift=S_t*M;
S_new = S_t+vola+vola'+(drift+drift'+drift_fix)*h;
S=cat(3,S,S_new);
eigS=[eigS;sort(eig(S_new)')];
timestep=[timestep;T];
end

When that is done, this simulates the **a single Euler-Maruyama realization of a path**:

function [X]=EulerSinglePath(T,h,S_0)
% simulation of Wishart based process
% parameters are drawn from Gauthier and Possamai paper
% Efficient Simulation
% dimension of the Wishart variance process
d = 2;
alpha =4.8;
% initial log spot
X_0 = log(S_0);
% Variance initial level
s_0 = [0.5, 0.; 0., 0.5];
% drift correction
Q = [0.35, 0.; 0., 0.4];
% drift component
M = [-0.6, 0.; 0., -0.4];
% correlation
R = [-0.5, 0.; 0.,-0.5];
% Simulation of wishart variance
[S,eigS,timestep,minh] = Wishart(T,h,d,alpha,Q,M,s_0);
% after simulation of matrix valued process, simulation of asset is given here
% precomputing the correlation structure
IminRR = eye(d,d)-R*R';
[V_ImRR,D_ImRR]=eig(IminRR);
sqImR=V_ImRR*sqrt(D_ImRR)*V_ImRR';
% ciclyng in the third dimension of the matrix
n_steps = size(S)(3);
X = zeros(1, n_steps);
% setting the initial log asset level
X(1) = X_0;
for n =2:n_steps
TrS_t = trace(S(:,:,n));
drift = -0.5*TrS_t*h;
% VOLATILITY
[V_t,D_t]=eig(S(:,:,n));
% the formula below corresponds to
% Z_t := W_tR^T + B_t \sqrt{I_d − RR^T}
W_t=sqrt(h)*normrnd(0,1,d,d);
B_t=sqrt(h)*normrnd(0,1,d,d);
Z = W_t * R' + B_t * sqImR;
vola=V_t*sqrt(D_t)*V_t'*Z;
X_old = X(n-1);
X(n) = X_old + drift + trace(vola);
X_old=X(n);
end

### NOTES

Rogers, Tehranchi, Can the implied volatility surface move by parallel shifts?