%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Week 9 Homework – Solution
% Part 1: Monte Carlo OLS (montecarlo_ols)
% Part 2: Numerical Integration of CRRA Utility
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear; close all; clc;

%% ------------------------------------------------------------------------
% PART 1 – Monte Carlo OLS
% -------------------------------------------------------------------------

% Parameters
beta0_true = 1;
beta1_true = 0.5;
sigma_true = 1;

R = 2000;                 % number of replications
n_values = [20, 50, 200]; % sample sizes to study

rng(123); % reproducibility

% Storage for results
means_beta1 = zeros(length(n_values),1);
vars_beta1  = zeros(length(n_values),1);

for j = 1:length(n_values)

    n = n_values(j);

    % Call the Monte Carlo function
    beta1_hat = montecarlo_ols(n, R, beta0_true, beta1_true, sigma_true);

    % Empirical mean and variance
    means_beta1(j) = mean(beta1_hat);
    vars_beta1(j)  = var(beta1_hat);

    % Display summary
    fprintf('--- Monte Carlo for n = %d ---\n', n);
    fprintf('Empirical mean of beta1_hat:     %.4f\n', means_beta1(j));
    fprintf('Empirical variance of beta1_hat: %.4f\n\n', vars_beta1(j));

    % Plot histogram
    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 approximation with empirical mean/variance
    mu_hat = means_beta1(j);
    se_hat = sqrt(vars_beta1(j));
    grid_min = mu_hat - 4*se_hat;
    grid_max = mu_hat + 4*se_hat;
    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

% (Optional) Compare empirical variances with 1/n
fprintf('Theoretical variance (approx) = 1/n (since Var(x)=1, sigma^2=1)\n');
for j = 1:length(n_values)
    fprintf('n = %d: empirical var = %.4f,  1/n = %.4f\n', ...
        n_values(j), vars_beta1(j), 1/n_values(j));
end
fprintf('\n');


%% ------------------------------------------------------------------------
% PART 2 – Numerical Integration of CRRA Utility
% -------------------------------------------------------------------------

% CRRA parameters and range
gamma = 2;
u = @(c) (c.^(1-gamma) - 1) ./ (1-gamma);  % gamma != 1

c_min = 0.1;
c_max = 10;

% Different grid sizes
N_values = [50, 200, 1000];

fprintf('--- Numerical Integration of CRRA Utility ---\n');

results_trap = zeros(length(N_values),1);
results_simp = zeros(length(N_values),1);

for j = 1:length(N_values)

    N = N_values(j);
    c_grid = linspace(c_min, c_max, N);
    uc = u(c_grid);

    % Trapezoidal rule
    int_trap = trapz(c_grid, uc);

    % Simpson's rule (manual implementation)
    h = c_grid(2) - c_grid(1);
    % Need an odd number of points for standard Simpson's rule:
    if mod(N,2) == 0
        % If N is even, drop the last point for Simpson
        uc_s = uc(1:end-1);
        h_s  = c_grid(2) - c_grid(1);
        N_s  = length(uc_s);
        int_simp = (h_s/3) * ( uc_s(1) + 4*sum(uc_s(2:2:end-1)) + ...
                               2*sum(uc_s(3:2:end-2)) + uc_s(end) );
    else
        % Direct Simpson if N is odd
        int_simp = (h/3) * ( uc(1) + 4*sum(uc(2:2:end-1)) + ...
                             2*sum(uc(3:2:end-2)) + uc(end) );
    end

    results_trap(j) = int_trap;
    results_simp(j) = int_simp;

    fprintf('N = %4d:  trapz = %10.6f   Simpson = %10.6f\n', ...
        N, int_trap, int_simp);
end

% Simple plot to illustrate convergence as N increases
figure;
plot(N_values, results_trap, '-o', 'LineWidth', 1.5); hold on;
plot(N_values, results_simp, '-s', 'LineWidth', 1.5);
xlabel('Number of grid points N');
ylabel('Approximate integral of u(c)');
title('Convergence of Numerical Integration (CRRA Utility)');
legend('Trapezoidal', 'Simpson', 'Location', 'best');
grid on;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Local function: montecarlo_ols
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function beta1_hat = montecarlo_ols(n, R, beta0, beta1, sigma)
% montecarlo_ols
%   Simulate the sampling distribution of the OLS slope in a simple
%   linear regression model:
%       y_i = beta0 + beta1 * x_i + u_i
%   where x_i ~ N(0,1), u_i ~ N(0, sigma^2).
%
%   Inputs:
%       n      - sample size
%       R      - number of Monte Carlo replications
%       beta0  - true intercept
%       beta1  - true slope
%       sigma  - standard deviation of the error term
%
%   Output:
%       beta1_hat - R-by-1 vector of estimated slopes

    beta1_hat = zeros(R,1);

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

        % Construct y
        y = beta0 + beta1 * x + u;

        % OLS regression y = b0 + b1 * x + e
        X = [ones(n,1), x];
        bhat = X \ y;
        beta1_hat(r) = bhat(2);
    end
end