%% Problem Sheet 7
% Methods of Monte Carlo Simulation

clear all; close all; clc;

%% Exercise 3 (importance sampling)

N = 10^5;
lambda = 1/2;
fprintf('True probability: %f\n\n', exp(-10*lambda));

% Standard Monte Carlo
Y_cmc = zeros(N,1);
for i = 1:N
    X = -(1/lambda)*log(1-rand);
    Y_cmc(i) = (X > 10);
end
fprintf('--- Standard Monte Carlo ---\n');
fprintf('The estimated probability is %f\n', mean(Y_cmc));
fprintf('The estimated relative error is %f\n\n', std(Y_cmc)/(sqrt(N)*mean(Y_cmc)));


% Variance Minimization Method (calculated theta)
theta_vm = (-1 + sqrt(1 + 100*lambda^2))/10;
f_div_g = @(x, theta) (lambda/(lambda-theta))*exp(-theta*x);

Y_vm = zeros(N,1);
for i = 1:N
    X = -(1/(lambda-theta_vm))*log(1-rand);
    Y_vm(i) = (X > 10)*f_div_g(X, theta_vm);
end
fprintf('--- Importance Sampling (calculated theta) ---\n');
fprintf('The optimal theta is %f\n', theta_vm);
fprintf('The estimated probability is %f\n', mean(Y_vm));
fprintf('The estimated relative error is %f\n\n', std(Y_vm)/(sqrt(N)*mean(Y_vm)));


% Variance Minimization Method (estimated theta)
N_init = 10^4;
X_init = -log(rand(N_init,1))/lambda;
X_sample = X_init(X_init > 10);
theta_est = fminsearch(@(theta) second_moment(theta, lambda, X_sample), 0);

X = -log(rand(N,1))/(lambda-theta_est);
Y_vm_est = (X > 10).*f_div_g(X, theta_est);

fprintf('--- Importance Sampling (estimated theta) ---\n');
fprintf('The estimated optimal theta is %f\n', theta_est);
fprintf('The estimated probability is %f\n', mean(Y_vm_est));
fprintf('The estimated relative error is %f\n\n\n', std(Y_vm_est)/(sqrt(N)*mean(Y_vm_est)));


%% Exercise 6 (incompetent business man)

mu = -30;
sigma = sqrt(10000);

X0 = 10000;
gamma = 11000;
N = 10^4;

% a)

company_sold = zeros(N,1);
for i=1:N
    % simulate until sold or bankrupt
    X = X0;
    while (X >= 0) && (X < gamma)
        Y = mu + sigma*randn;
        X = X + Y;
    end
    % check why we quit the loop (sold or bankrupt)
    company_sold(i) = (X >= gamma);
end

fprintf('--- Standard Monte Carlo ---\n');
fprintf('The estimated probability is %f\n', mean(company_sold));
fprintf('The estimated standard deviation of the estimator is %e\n\n', std(company_sold)/sqrt(N));

% b)

mu_is = -mu;

company_sold_is = zeros(N,1);
for i=1:N
    % simulate until sold or bankrupt
    X = X0;
    f_div_g = 1;
    while (X >= 0) && (X < gamma)
        Y = mu_is + sigma*randn;
        X = X + Y;
        f_div_g = f_div_g * normpdf(Y, mu, sigma) / normpdf(Y, mu_is, sigma);
    end
    company_sold_is(i) = (X >= gamma)*f_div_g;
end

fprintf('--- Importance Sampling (calculated theta) ---\n');
fprintf('The estimated probability is %f\n', mean(company_sold_is));
fprintf('The estimated standard deviation of the estimator is %e\n\n', std(company_sold_is)/sqrt(N));