Replicate a density plot?! How To?

Hey there,

I have read the paper "The Heston Model A Practical Approach von Moodley (2005)". On page 3 there is the following plot of a density. It shows how the density changes, with a change in the parameter rho:

The Code given is the following:

S0 = 10; V0 = .01; r = 0; k = 2; 
theta =.01; sigma= .1; delT = .02; rho =-0.9; 
numberOfSimulations = 1000000; 

i = 1:numberOfSimulations; 
NormRand1 = randn(1,numberOfSimulations); 
NormRand2 = randn(1,numberOfSimulations); 
S = zeros(1,numberOfSimulations); 
V = zeros(1,numberOfSimulations); 
V(i) = V0 + k*(theta - V0)*delT + sigma*sqrt(V0)* ... 
(rho*NormRand1 + sqrt(1- rho^2)*NormRand2)*sqrt(delT); 
V = abs(V); %prevents negative volatilities 
S(i) = S0 + r*S0*delT + S0*V.^(0.5).*NormRand1*sqrt(delT);

Nothing more is given than this code. But i dont get how to build it, with this code? Maybe someone can help me with it.

kind regards david

Answer 1

[Geoff Hayes]

David - remember what the density function is. If you look at the minimum and maximum of the S array you will notice that these values (respectively) fall around 9.5 and 10.5. What the author probably did was to create a histogram of the data (from S) and then plot the results to get the figure. Something like

 figure;
 hold on;

 S0 = 10; V0 = .01; r = 0; k = 2; 
 theta =.01; sigma= .1; delT = .02; 
 numberOfSimulations = 1000000; 
 i = 1:numberOfSimulations; 
 NormRand1 = randn(1,numberOfSimulations); 
 NormRand2 = randn(1,numberOfSimulations); 
 S = zeros(1,numberOfSimulations); 
 V = zeros(1,numberOfSimulations); 

 for k=1:3
    switch k
        case 1
            rho = 0.9;
            ls  = '-k';
        case 2
            rho = 0.0;
            ls  = '--r';
        case 3
            rho = -0.9;
            ls  = '-.b';
    end

    V(i) = V0 + k*(theta - V0)*delT + sigma*sqrt(V0)* ... 
    (rho*NormRand1 + sqrt(1- rho^2)*NormRand2)*sqrt(delT); 
    V = abs(V); %prevents negative volatilities 
    S(i) = S0 + r*S0*delT + S0*V.^(0.5).*NormRand1*sqrt(delT);

     [nelements,centers] = hist(S,50);

     plot(centers,nelements/sum(nelements),ls);
 end

The majority of the above code is what you already have. The for loop allows us to calculate S for different rho and then plot the results. Note the

 [nelements,centers] = hist(S,50);

sorts the data in 50 bins whose centres are given by centers with nelements per bin. We then plot the result as

 plot(centers,nelements/sum(nelements),ls);

The above divides by sum(elements) to give a reasonable approximation to the figure.

[Geoff Hayes]

David's answer moved here

Hi Geoff,

I can't test it now but after work I will. Thank you so so much in advance :)! You See I am a newbie with Matlab (Never done it before). But with a little help I will get through it.

Regards David