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

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