How to integrate a shifted lognormal distributed random variable

Dear all,

I'm quite new to Matlab and struggeling to integrate a continuous random variable by two parts. In the beginning I assume W=1+w is lognormally distributed with mean = 1 and standard deviation of 0.05. The mean is specifically chosen such that w has a zero mean and has a support of [-1, Inf). I tried to integrate a hard coded version of a lognormal density. This doesn't give me a useful result since log is not defined for -1 yet I tried the same procedure with W the result still doesn't make sense.

% Generating random variable W

W = lognrnd(0, 0.05, 1, 2000);

% Shifting W to have lognormal random variable w with zero mean and
% has a support of [-1, Inf)

w = W - 1;

% Transforming omega to s_w = a*w with a = 9.0083

s_w = 9.0083*w;

% My aim is now to integrate s_w from -1 to 0 and from 0 to Inf
% Since I can only integrate a continous random variable numerically I
% tried to hard code the density and use the integral function

%mu = mean(s_w);
%sigmasq = var(s_w);

%s_w = @(x) exp(-(log(x) - mu).^2./(2.*sigmasq)./(x.*sqrt(2.*pi.*sigmasq)));

%S_neg = integral(@(x)s_w(x), -1, 0);
%S_pos = integral(@(x)s_w(x), 0, Inf);

%d = -S_neg/S_pos;

% This does obviously not work because of the support of log yet I actually
% need to integrate

% Even my approach of tranforming W to s(W) doesn't work since both values
% are close to zero which they shouldn't be according to a density plot

s_W = 9.0083*W;

mu = mean(s_W);
sigmasq = var(s_W);

s_W = @(x) exp(-(log(x) - mu).^2./(2.*sigmasq)./(x.*sqrt(2.*pi.*sigmasq)));

S_neg = integral(@(x)s_W(x), (-1+mu), mu);
S_pos = integral(@(x)s_W(x), mu, Inf);

d = -S_neg/S_pos;

I hope I gave all the necessary information.

Answer 1

[Matt J]

For the edited version of your post,

a=-1+mu;
b=mu;
sigma=sqrt(sigmasq);

Sneg = logncdf(b,mu,sigma) - logncdf(a,mu,sigma);
Spos = logncdf(b,mu,sigma,'upper');

[Lukas Huppertz]

Thank you but why is Sneg still zero and Spos = 1, when I take a look at the density plot of s_W this result seems impossible to me.

[k,pdf_s_W] = ksdensity(s_W)
plot(pdf_s_W,k)

[Matt J]

Because you are still stuck on the idea that

@(x) exp(-(log(x) - mu).^2./(2.*sigmasq)./(x.*sqrt(2.*pi.*sigmasq)));

is the pdf of 9.0083*W. It is not.

If the goal is to determine the probabilities of

Low<=9.0083*W <=High

    9.0083*W>=High

then that is the equivalent to

    Low/9.0083 <= W <=High/9.0083

        W>=High/9.0083

Because W is lognormal with parameters mu=0 and sigma=.05, this is readily adaptable from what I've already proposed above:

a=Low/9.0083;
b=High/9.0083;
mu=0;
sigma=.05;

Sneg = logncdf(b,mu,sigma) - logncdf(a,mu,sigma);
Spos = logncdf(b,mu,sigma,'upper');

[Lukas Huppertz]

If

@(x) exp(-(log(x) - mu).^2./(2.*sigmasq)./(x.*sqrt(2.*pi.*sigmasq)));

isn't the pdf of 9.0083*W what is the pdf of it? My goal is the probabilities or areas beneath the density above. And since my low border is always in the same area (in the graph above it would probably be 8) and the high border is correspondingly about +1 higher, I cannot see why Sneg is 0 and Spos is 1 instead of values like Sneg is 0.35 and Spos is 0.6.

[Matt J]

Well, since you've randomly generated mu and sigma, I can't reproduce what you see precisely. However, when I run your code, I obtain

    mu =

      9.0209

   sigma =

      0.4427

When I plug these parameters, as you have done, into the formula for a lognormal pdf and then plot it,

s_W = @(x) exp(-(log(x) - mu).^2./(2.*sigmasq))./( x.*sqrt(2.*pi.*sigmasq) );

   fplot(s_W,[0,15000]);

the curve below is obtained. So as you can see, the distribution has essentially no probability mass for x<mu. Thus, you get Sneg=0.

This will always be the case for mu significantly greater than 1. Remember that the mean of a lognormal is exp(mu+sigma^2/2), which grows very fast with mu.

If ...isn't the pdf of 9.0083*W what is the pdf of it?

Since W is Lognormal(0,0.05) then c*W will be Lognormal(log(c), .05).

[Lukas Huppertz]

Wow thank you that helped alot. So, if my approach to make 9.0083*w or 9.0083*W integratable with coding the lognormal density as a function does not work, how can I integrate 9.0083*w or 9.0083*W? In all my replications the shape is quite similar to the one above and that is basically my integration goal.

And doesn't this linear transformation also influence the variance besides the mean?

[Matt J]

Well, 9.0083*W is just another Lognormal with mu=log(9.0083). You already know how to integrate those using logncdf. w is a shifted verson of this, but instead of shifting the distribution, you can just shift the limits of integration.

[Lukas Huppertz]

Thank you Matt. Even though I can't see how it solves my problem with w and its positive probabilities for negative values, you helped me alot figuring out how to compute the integral itself which is the first, more important step.

Answer 2

[Matt J]

Maybe you can explain the ulterior motive of this. The s_om that you've defined has the form of a non-shifted lognormal pdf, but w_s is a shifted/scaled lognormal. Are you sure you shouldn't be integrating the corresponding shifted/scaled pdf?

Incidentally also, the shifting w = W - 1 does not make w zero-mean. The mean of W as you've constructed it is exp(0+.05^2/2), so the mean of w will be exp(0+.05^2/2)-1.

[Lukas Huppertz]

Thank you for pointing out the mistakes in my code, I've edited the code to make it more clear. Basically, I want to integrate s_w but I simply don't know how to. That is why I tried like this.

Also, thank you for pointing out the mistake about the mean, I will try to do it correctly. My problem is that if I calculate the needed mean it should be -0.00125 but this will not fit for only 2000 random draws. I use 2000 draws since I'm trying to replicate something and the authors use 2000 trapzoids for a numerical integration and as far as I've found out I can only draw 2000 times such that both vector lengths are the same.

The motive is that s_w is the distribution of banks with a excess or deficit based on the shock w and I need the relative proportion of banks with excess to banks with deficit

[Matt J]

Well, either you or the paper you're following has to have made a mathematical mistake in the problem formulation. As you've posed the problem now, the result is trivial:

 S_neg = 0;
 S_pos = 1;

This is because s_W is the pdf of a (non-shifted) lognormal distribution, so it's integral from 0 to Inf has to be 1 (lognormal variables are positive with probability 1). Likewise, any integral outside of this range has to be 0.

[Lukas Huppertz]

Thank you Matt. I've made the mistake that I didn't shift the boundaries from s_w and to the corresponding integral boundaries for s_W. I also noticed that the result for d gets more realistic when I change the upper boundary of S_pos from Inf to a more realistic number like 11. This does make sense since the integral function may have trouble to specify the step length of the numeric approximation steps. I simply do not understand why these values for S_neg are still so low because when I compute it on my own the value is way higher.