simple integration of joint normal - TPT doesn't work

Pete sherer — Tue, 17 Mar 2009 18:04:01 -0400

Hi,

This problem seems to be simple and should have any problem when coding. However, it seems to create problems on my side. Below, I am trying to obtain a marginal distribution from a joint lognormal distribution. I howver can't get the same answer as the theoretical one when correlation is different from zero.

levSde2Exc = 1;
trho = 0.0;
[ meanSde1, meanSde2 ] = deal( log(1), log(2));
[ var1, var2] = deal( (0.6), (0.8));

% define integration range over the first variable (i.e., Sde_1)
min1 = meanSde1 - 5*sqrt(var1); % min( log(Sde_1levIM));
max1 = meanSde1 + 5*sqrt(var1); % max( log(Sde_1levIM));

Sde_1levIM = exp( linspace( min1, max1, 3000) );
numSdeT1 = length( Sde_1levIM);
fsde2_condSde1 = zeros( 1, numSdeT1);

% using the Total probability Theorem to integrate jointly lognormal distribution to obtain the marginal of the 2nd variable (i.e., Sde_2)
fxdx = normpdf( log( Sde_1levIM), meanSde1, sqrt( var1));
fxdx = fxdx./ sum( fxdx);
for runSde1 = 1: numSdeT1 % I can collapse this for loop
    levSde1 = Sde_1levIM( 1, runSde1);
    mean_22Cond = meanSde2 + sqrt( var2)/sqrt(var1)*trho*( levSde1 - meanSde1);
    var_22Cond = (1-trho^2)* var2;

    tpSde2_condX_Sde1 = normpdf( log( levSde2Exc), mean_22Cond, sqrt( var_22Cond));
    fsde2_condSde1( 1, runSde1) = tpSde2_condX_Sde1.* fxdx( 1, runSde1);
end % for runSde1 = 1: length(scalarPSHA.Sde_1.levIM)

%% Integrated fsde2
fsde2_cal = sum( fsde2_condSde1, 2);
%% Marginal fsde2
fsde2_theo = normpdf( log( levSde2Exc), meanSde2, sqrt( var2) );

when correlation (trho = 0), I get the same answer between fsde2_cal and fsde2_theo which is 0.3303.

when trho = 0.5, fsde2_cal = 0.1342 which is different from 0.3303.

Any suggestions to this simple problem. Thank you so much in advance.

Re: simple integration of joint normal - TPT doesn't work

Pete sherer — Tue, 17 Mar 2009 18:20:17 -0400

found the BUG!!
I forgot to take log of levSde1 when calculating mean_22Cond.

it should read
mean_22Cond = meanSde2 + sqrt( var2)/sqrt(var1)*trho*( log(levSde1) - meanSde1)

MATLAB Central Newsreader - simple integration of joint normal - TPT doesn't work

simple integration of joint normal - TPT doesn't work

Re: simple integration of joint normal - TPT doesn't work