%% Problem Sheet 6
% Methods of Monte Carlo Simulation

clear all; close all; clc;

%% Exercise 2 (estimation of pi)

N = 10^4;

% a)

U_1 = rand(N,1);
X_a = 4*sqrt(1 - U_1.^2);

fprintf('The estimate of pi using conditional MC is %f.\n', mean(X_a));

% b)

U = rand(N/2, 2);
Y = 4*(U(:,1).^2 + U(:,2).^2 <= 1);
Y_star = 4*((1-U(:,1)).^2 + (1-U(:,2)).^2 <= 1);
X_b = [Y; Y_star];

fprintf('The estimate of pi using antithetic variables is %f.\n', mean(X_b));

% Justification of the estimator:
% 
% If (U1, U2) is in the quarter of the circle, (1-U1, 1-U2) is likely to be
% not. The other way round, if (U1, U2) is not in the circle, (1-U1, 1-U2)
% surely is.
%
% => The two indicators have to be negatively correlated, the constant does
% not change this.

% c)

fprintf('The estimated variances of the estimators are %e for a) and %e for b).\n\n', ...
        var(X_a)/N, var(X_b)/N);
% => Conditional MC performs better in this case.


%% Exercise 4 (random walk)

close all; clear all;

N = 10^4;

% a) (crude Monte Carlo)

% generate Z_i
U = rand(N, 10);
Z = zeros(N,10);
Z(U < 0.5) = 1;
Z(U >= 0.5) = -1;

% check if sums are > 1
Y = (sum(Z, 2)>1);
fprintf('Estimated probability using standard MC: %f\n', mean(Y));

% b)

% generate Z_i for N/2 random walks
U2 = rand(N/2, 10);
Z2 = zeros(N/2, 10);
Z2(U2 < 0.5) = 1;
Z2(U2 >= 0.5) = -1;

% check if X_10 (the sum of Z_1, ..., Z_10) is > 1 ... 
Y2 = (sum(Z2, 2) > 1); % ... for the original random walks
Y2_star = (-sum(Z2, 2) > 1); % ... and for their reflections
fprintf('Estimated probability using antithetic variables: %f\n\n', mean([Y2; Y2_star]));

% c) (Comparison)

var_cmc = var(Y)/N;
var_av = var(Y2 + Y2_star)/(2*N);
fprintf('Variance of the estimator using standard MC: %d\n', var_cmc);
fprintf('Variance of the estimator using antithetic variables: %d\n', var_av);
fprintf('Quotient of the two: %f\n\n', var_av/var_cmc);


%% Exercise 6 (stratified sampling)

strata = [0, 1; 1, 2; 2, 6; 6, 10];
f = @(x) (9/50)*x.^3 -(7/2)*x.^2 + 20*x -22;

% a)

N = 10^3;
U = rand(N,1)*10;

p = zeros(size(strata,1), 1);
sigma = zeros(size(strata,1), 1);

for i=1:size(strata,1)
    X = U((U >= strata(i,1)) & (U < strata(i,2)));
    p(i) = length(X) / N;
    sigma(i) = std(10*f(X));
end
    
q = sigma.*p / sum(sigma.*p);

fprintf('The estimated proportions are: %f, %f, %f, %f\n', q);
fprintf(['With a total sample size of 10^5, this corresponds to sample', ...
        ' sizes of: %d, %d, %d, %d\n\n'], round(q*10^5));

% b)

N = 10^5;
N_i = round(q .* N);

Y_bar = zeros(size(strata,1),1);
sigma_sq = zeros(size(strata,1),1);

for i=1:size(strata,1)
    U = rand(N_i(i),1)*(strata(i,2)-strata(i,1)) + strata(i,1);
    Y_bar(i) = mean(10*f(U));
    sigma_sq(i) = var(10*f(U))/N_i(i);
end
    
est_int = sum(p.*Y_bar);
est_var = sum((p.^2).*sigma_sq);

fprintf('---Stratified Sampling---\n');
fprintf('The estimated value of the integral is %f.\n', est_int);
fprintf('The estimated variance is %f.\n\n', est_var);
    
% c)

U = rand(N,1)*10;
est_int_2 = mean(10*f(U));
est_var_2 = var(10*f(U))/N;

fprintf('---Standard Monte Carlo---\n');
fprintf('The estimated value of the integral is %f.\n', est_int_2);
fprintf('The estimated variance is %f.\n\n', est_var_2);