%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Week 9 In-Class Challenge – Solution File (Advanced Monte Carlo)
% Numerical Differentiation, Integration, and Monte Carlo for OLS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear; close all; clc;

%% ------------------------------------------------------------------------
% PART 1 – Numerical Differentiation and Integration
% -------------------------------------------------------------------------

%% Task A: Numerical derivative of f(x) = log(x) at x = 2

f  = @(x) log(x);
x0 = 2;

h = 1e-5;    % step size

% Analytical derivative
fp_true = 1 / x0;

% Forward difference
fp_forward = (f(x0 + h) - f(x0)) / h;

% Central difference
fp_central = (f(x0 + h) - f(x0 - h)) / (2*h);

disp('--- Task A: Numerical Differentiation ---');
fprintf('Analytical derivative:      %.6f\n', fp_true);
fprintf('Forward difference approx:  %.6f\n', fp_forward);
fprintf('Central difference approx:  %.6f\n\n', fp_central);


%% Task B: Numerical integration of ∫_1^3 x^2 dx

g = @(x) x.^2;

a = 1;
b = 3;
N = 200;                       % number of grid points
x = linspace(a, b, N);
fx = g(x);

% Analytical integral
int_true = (b^3 - a^3) / 3;

% Trapezoidal rule
int_trap = trapz(x, fx);

% Simpson's rule
h = x(2) - x(1);
int_simp = (h/3) * ( fx(1) + 4*sum(fx(2:2:end-1)) + ...
                     2*sum(fx(3:2:end-2)) + fx(end) );

disp('--- Task B: Numerical Integration ---');
fprintf('Analytical integral:        %.6f\n', int_true);
fprintf('Trapezoidal approximation:  %.6f\n', int_trap);
fprintf('Simpson approximation:      %.6f\n\n', int_simp);


%% ------------------------------------------------------------------------
% PART 2 – Advanced Monte Carlo: OLS Slope Distribution and Coverage
% -------------------------------------------------------------------------

% Data-generating process:
% y_i = beta0 + beta1 * x_i + u_i
beta0_true = 1;
beta1_true = 0.5;

R  = 2000;        % number of replications
n_vec = [50, 200]; % two sample sizes

rng(123);        % for reproducibility

% Pre-allocate storage
beta1_hat_all = cell(length(n_vec), 1);
cov_all       = cell(length(n_vec), 1);   % coverage indicators
mean_emp      = zeros(length(n_vec), 1);
var_emp       = zeros(length(n_vec), 1);
cov_rate      = zeros(length(n_vec), 1);

for j = 1:length(n_vec)

    n = n_vec(j);
    beta1_hat = zeros(R,1);
    coverage  = false(R,1);

    for r = 1:R
        % Draw regressors and errors
        x = randn(n,1);          % N(0,1)
        u = randn(n,1);          % N(0,1)

        % Construct y
        y = beta0_true + beta1_true * x + u;

        % OLS regression: y = b0 + b1 * x + error
        X = [ones(n,1), x];      % design matrix
        bhat = X \ y;            % OLS via backslash

        beta1_hat(r) = bhat(2);

        % Compute usual homoskedastic standard error of b1
        yhat = X * bhat;
        residuals = y - yhat;
        s2 = sum(residuals.^2) / (n - 2);                 % sigma^2 hat
        Sxx = sum( (x - mean(x)).^2 );
        se_b1 = sqrt( s2 / Sxx );

        % 95% CI for beta1
        lower = beta1_hat(r) - 1.96 * se_b1;
        upper = beta1_hat(r) + 1.96 * se_b1;

        coverage(r) = (beta1_true >= lower) && (beta1_true <= upper);
    end

    % Store results
    beta1_hat_all{j} = beta1_hat;
    cov_all{j}       = coverage;

    % Empirical statistics
    mean_emp(j) = mean(beta1_hat);
    var_emp(j)  = var(beta1_hat);
    cov_rate(j) = mean(coverage);  % share of CIs containing true beta1
end

% Display summary statistics
disp('--- Part 2: Monte Carlo for OLS Slope ---');
for j = 1:length(n_vec)
    n = n_vec(j);
    fprintf('Sample size n = %d\n', n);
    fprintf('  Empirical mean of beta1_hat:     %.4f\n', mean_emp(j));
    fprintf('  Empirical variance of beta1_hat: %.4f\n', var_emp(j));
    fprintf('  Empirical 95%% CI coverage rate:  %.4f\n\n', cov_rate(j));
end

% Theoretical variance for beta1_hat under homoskedasticity:
% Var(beta1_hat) ≈ sigma^2 / (n * Var(x)), here sigma^2 = 1, Var(x) = 1
fprintf('Theoretical variance (approx) = 1/n\n');
fprintf('For n = 50:  1/n = %.4f\n', 1/50);
fprintf('For n = 200: 1/n = %.4f\n\n', 1/200);

%% Plots: histograms for each sample size

for j = 1:length(n_vec)
    n = n_vec(j);
    beta1_hat = beta1_hat_all{j};

    figure;
    histogram(beta1_hat, 40, 'Normalization', 'pdf'); hold on;
    title(sprintf('Sampling Distribution of \\beta_1 (n = %d)', n));
    xlabel('\beta_1^{hat}');
    ylabel('Density');

    % Overlay normal density with empirical mean and variance
    mu_hat = mean(beta1_hat);
    se_hat = sqrt(var(beta1_hat));

    grid_min = mu_hat - 4*se_hat;
    grid_max = mu_hat + 4*se_hat;
    if grid_min < min(beta1_hat), grid_min = min(beta1_hat); end
    if grid_max > max(beta1_hat), grid_max = max(beta1_hat); end

    bgrid = linspace(grid_min, grid_max, 400);
    plot(bgrid, normpdf(bgrid, mu_hat, se_hat), 'LineWidth', 2);

    legend('Monte Carlo histogram','Normal approx','Location','best');
end