I'm converting some calculations related to an earlier question to MATLAB.
P LB UB ExpLoss
Event 1 2.0% $40,000 $200,000 $2,016
Event 2 5.0% $500,000 $2,000,000 $54,643
Event 3 10.0% $400,000 $2,500,000 $116,785
Event 4 15.0% $100,000 $5,000,000 $215,104
Four Events have the estimated probability, a financial consequence (measured by estimated lower and upper bounds) and an expected loss - based on the assumption of a 90% "confidence interval" and an assumed lognormal distribution.
The expected value, as r.e.s. explained, is calculated by
$ \mathbb{E}[\text{ loss }]=P(E)\,\mathbb{E}[X] = P(E)\,\exp\left({\mu+\frac{1}{2}\sigma^2}\right) $
I have constructed some MATLAB code intended to capture this relationship for the four events. In it, I compute the expected Loss, and get the same results that Hubbard does in his spreadsheet (see previous question if needed). So far so good. However, I then tried to recreate that expectation via simulation. For a large number of trials, I plug a random number into the calculation and compute the expected loss for that trial. I'm assuming that the value of "expected loss" means "the amount that would be lost, on average for a large number trials". But the output of my code is as follows:
Expected Loss fprintf ('%.2f\n', Event) fprintf('%.2f\n', accumulator/numruns)
2016.27 1006.87 20.01
54642.57 27298.25 1350.85
116784.55 58527.96 5814.29
215104.47 107859.02 16115.65
The "if" statement conditional is executing the correct number of times for each of the four events (approx 2%, 5%, 10%, and 15%)
I thought that dividing the accumulated loss by the number of trials (numruns) should approximate the expected value. To my surprise, the accumulated loss divided by the much smaller non-zero events (counter, as seen in the second column) is approximating half of the expected value.
So, my question is this: Is my understanding of the Expected Loss incorrect? or is my code incorrect? (or both?!). Please let me know if I have failed to explain everything clearly enough, and thanks in advance.
%%% demo of risk
LB = [40000 500000 400000 100000];% These are the lower and upper bounds
UB = [200000 2000000 2500000 5000000]; %for the four events
P = [.02 .05 .1 .15];% The probabilities of the four events
numruns = 1000000; %the number of iterations we intend to do
logUB = log(UB); % The natural logs of the upper and lower bounds
logLB = log(LB);
mu = (logUB + logLB)/2; %the mean and variance of the bounds. The 90 percent
sigma = (logUB-logLB)/3.28971;%confidence interval is established by 3.2 stds
expLoss = P.*exp((mu + sigma.^2/2)); % This confirms that we're getting
%the same value as the spreadsheet
fprintf('%8.2f\n',expLoss);
randP = rand(1,numruns);%this creates an array of random probabilities
Event = zeros(1,4); %an array for the four events
counter = zeros(1,4); %This counts the # times the "if" statement executes
accumulator = zeros(1,4); %This gathers the loss incurred from each of the four events
for i = 1:4 % We do one run of the inner loop for each of the four events
for k =1:numruns
==> if randP(k)< P(i) %If the Event occurs, compute the cost and add it
accumulator(i) = accumulator(i) + randP(k)*(exp((mu(i) + sigma(i).^2/2)));
counter(i) = counter(i) + 1;
end
end
Event(i)= accumulator(i)/counter(i); %This is the average loss per "non-zero" event.
end
fprintf('\n\n');
fprintf('%.2f\n', Event);
fprintf('\n\n');
fprintf('%.2f\n',accumulator/numruns)