Market completion with Wishart variance

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:

S_0 1.
T 1.
\Delta t 0.01.
d 2
\alpha 4.8

wishart_paths_20

 

wishart_paths_100

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?

Octave meshgrid, surf and video slideshow

Since the release of the splendid GUI, Octave has become again one of my favorite tools. Here is a simple and hopefully useful application. The following code is quite self-explanatory:

tx = ty = linspace (-4, 4, 41)';
[xx, yy] = meshgrid (tx, ty);
for ii=1:23
  alpha = 0.01*ii;
  tz = exp(-alpha*(xx .^2 + yy .^2 ));
  surf(tx, ty, tz);
  tit = strcat ("alpha=", num2str(alpha));
  title(tit);
  fname = sprintf("anim_%02i.jpg", ii);
  print (fname)
  ans = yes_or_no ("prompt")
endfor

 

After that, having installed ffmpeg, just issue from either bash or cmd shell:

ffmpeg.exe  -framerate 2 -start_number 1 -i "anim_%02d.jpg" -vcodec mpeg4 evolution.mp4

 

to obtain: