%% week6_challenge_starter.m
% Week 6 — In-class Challenge: Calibrate Solow savings rate s
% Goal: Find s that matches an observed output path y_data.
% Steps: (i) manual grid search; (ii) (bonus) fminsearch with reparam.

clear; close all; clc;

% Ensure Figures folder exists
if ~exist('Figures','dir'), mkdir('Figures'); end

%% ------------------ Model & Data Setup ------------------
% Parameters (fixed)
params.alpha = 0.33;   % capital share
params.delta = 0.08;   % depreciation
params.n     = 0.01;   % population growth
params.g     = 0.02;   % technology growth
params.A0    = 1.0;    % initial technology
params.K0    = 1.0;    % initial capital
T            = 40;     % horizon

% Synthetic "observed" data generated with s_true + small noise
s_true   = 0.24;
y_model0 = solow_simulate(s_true, params, T);  % true model path
rng(123);
y_data   = y_model0 .* exp(0.05*randn(size(y_model0))); % ~noise-free for class

% Objective: SSE between model and data for a given s
obj = []; %% <-- Students to fill

%% ------------------ (1) Manual Grid Search ------------------
s_grid = linspace(0.05, 0.5, 200);
SSE    = zeros(size(s_grid));
for i = 1:length(s_grid)
    [] %% <-- Students to fill!
end
%% Students to find the minimium of SSE and record the corresponding index grid into s_hat_grid
s_hat_grid=[];   %%<-- Students to fill!

fprintf('Grid search estimate: s_hat = %.4f (true s = %.4f)\n', s_hat_grid, s_true);

% Plot objective over grid
[];   %%<-- Students to fill!

%% ------------------ (2) BONUS: fminsearch with reparam ------------------
% Re-parameterize s = sigmoid(x) to enforce 0 < s < 1
sigmoid = @(x) 1./(1 + exp(-x));
obj_x   = @(x) obj(sigmoid(x));

x0 = 0;  % corresponds to s0 = 0.5
[x_hat, fval] = []; %%<-- Students to fill!
s_hat_nelder  = sigmoid(x_hat);

fprintf('fminsearch estimate (reparam): s_hat = %.4f, SSE = %.4f\n', s_hat_nelder, fval);

% Compare fitted vs data
y_fit = solow_simulate(s_hat_nelder, params, T);
[];   %%<-- Students to fill!

%% ------------------ Helper: Solow simulator ------------------
function y = solow_simulate(s, P, T)
    % Clamp s to (0,1) for numerical safety in class
    s = max(min(s, 0.999), 0.001);
    K = zeros(T,1); A = zeros(T,1); Y = zeros(T,1);
    K(1) = P.K0; A(1) = P.A0;
    for t = 1:T
        Y(t) = K(t)^P.alpha * A(t)^(1-P.alpha);         % Y_t = K^alpha (A)^(1-alpha), L normalized
        if t < T
            K(t+1) = (1 - P.delta)*K(t) + s*Y(t);
            A(t+1) = (1 + P.g)*A(t);
        end
    end
    y = Y;  % observed output path (per effective worker scaling absorbed by A)
end