Path: news.mathworks.com!not-for-mail
From: <HIDDEN>
Newsgroups: comp.soft-sys.matlab
Subject: simple integration of joint normal - TPT doesn't work
Date: Tue, 17 Mar 2009 18:04:01 +0000 (UTC)
Organization: AIR Worldwide Corp
Lines: 39
Message-ID: <gpoomh$k8v$1@fred.mathworks.com>
Reply-To: <HIDDEN>
NNTP-Posting-Host: webapp-05-blr.mathworks.com
Content-Type: text/plain; charset="ISO-8859-1"
Content-Transfer-Encoding: 8bit
X-Trace: fred.mathworks.com 1237313041 20767 172.30.248.35 (17 Mar 2009 18:04:01 GMT)
X-Complaints-To: news@mathworks.com
NNTP-Posting-Date: Tue, 17 Mar 2009 18:04:01 +0000 (UTC)
X-Newsreader: MATLAB Central Newsreader 1076536
Xref: news.mathworks.com comp.soft-sys.matlab:525595


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.