%% week6_homework_starter.m
% Week 6 — Homework: Calibrate Solow savings rate s (Grid vs fminsearch vs fmincon)
% Replace the synthetic y_data if you have a provided series.

clear; close all; clc;
if ~exist('Figures','dir'), mkdir('Figures'); end

%% Parameters & (replaceable) data
P.alpha = 0.33; P.delta = 0.08; P.n = 0.01; P.g = 0.02; P.A0 = 1; P.K0 = 1;
T = 80;

% If you have an external series, load it here (uncomment and adapt):
% y_data = readmatrix('solow_output_data.csv'); y_data = y_data(:);

% For now, generate synthetic data with noise (students may replace):
s_true = 0.24;
rng(2025);
y_clean = solow_simulate(s_true, P, T);
y_data  = y_clean .* exp(0.05*randn(size(y_clean))); % small multiplicative noise

% Objective
obj = @(s) sum((solow_simulate(s, P, T) - y_data).^2);

%% (1) Manual grid search
tic;
s_grid = linspace(0.05, 0.5, 300);
SSE = arrayfun(obj, s_grid);
[~, i0] = min(SSE); s_hat_grid = s_grid(i0);
time_grid=toc;
%% (2) fminsearch with reparam s = sigmoid(x)
tic;
sigmoid = @(x) 1./(1 + exp(-x));   obj_x = @(x) obj(sigmoid(x));
[x_hat, fval1] = fminsearch(obj_x, 0);
s_hat_nelder = sigmoid(x_hat);
time_fminsearch=toc;

%% (3) fmincon with bounds 0 <= s <= 1
tic;
s0 = 0.3; lb = 0; ub = 1;
opts = optimoptions('fmincon','Display','off');
[s_hat_con, fval2] = fmincon(@(s) obj(s), s0, [], [], [], [], lb, ub, [], opts);
time_fmincon=toc;
% Display computation times for each method
fprintf('Time (Grid search): %.4f seconds\n', time_grid);
fprintf('Time (fminsearch): %.4f seconds\n', time_fminsearch);
fprintf('Time (fmincon): %.4f seconds\n', time_fmincon);
fprintf('Grid:      ŝ = %.4f\n', s_hat_grid);
fprintf('fminsearch ŝ = %.4f, SSE = %.4f\n', s_hat_nelder, fval1);
fprintf('fmincon    ŝ = %.4f, SSE = %.4f\n', s_hat_con,   fval2);

%% Plot objective and mark optima
fig1 = figure('Color','w','Position',[100 100 760 500]);
plot(s_grid, SSE, 'LineWidth',1.8); hold on;
ymin = min(SSE); ymax = max(SSE);
plot(s_hat_grid, interp1(s_grid,SSE,s_hat_grid), 'ko','MarkerSize',7,'LineWidth',1.5);
plot(s_hat_nelder, obj(s_hat_nelder), 'rs','MarkerSize',7,'LineWidth',1.5);
plot(s_hat_con, obj(s_hat_con), 'bd','MarkerSize',7,'LineWidth',1.5);
legend('Objective','Grid min','fminsearch','fmincon','Location','best');
xlabel('Savings rate s'); ylabel('SSE');
title('Objective vs s (Grid, fminsearch, fmincon)'); grid on; box on;
print(fig1, fullfile('Figures','fig_solow_obj_homework.pdf'), '-dpdf', '-r300');
print(fig1, fullfile('Figures','fig_solow_obj_homework.png'), '-dpng', '-r300');

%% Plot fit
y_fit = solow_simulate(s_hat_con, P, T);  % choose one (e.g., fmincon)
fig2 = figure('Color','w','Position',[100 100 760 500]);
plot(1:T, y_data, '-','LineWidth',1.8); hold on;
plot(1:T, y_fit,  '--','LineWidth',1.8);
xlabel('t'); ylabel('Output y_t'); title('Calibration Fit: Data vs Model');
legend('Data','Model (ŝ)','Location','best'); grid on; box on;
print(fig2, fullfile('Figures','fig_solow_fit_homework.pdf'), '-dpdf', '-r300');
print(fig2, fullfile('Figures','fig_solow_fit_homework.png'), '-dpng', '-r300');

%% Helper
function y = solow_simulate(s, P, T)
    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);
        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;
end