clear; close; clc;

alpha = 0.3;
k = 5;              % evaluation point
h = 1e-4;

f  = @(k) k.^alpha;
fp = alpha * k^(alpha-1);

fd = (f(k+h)-f(k))/h;
bd =  (f(k)-f(k-h))/h;
cd = (f(k+h)-f(k-h))/(2*h);

%%
clear; close; clc;
r = 0.04;
T = 30;
y = @(t) 10 + sin(t);     % example income profile

t = linspace(0,T,500);
f = y(t) .* exp(-r*t);
plot(t,f)

pv_trap = trapz(t,f);

% Simpson
h = t(2)-t(1);
pv_simp = (h/3) * (f(1) + 4*sum(f(2:2:end-1)) + ...
                   2*sum(f(3:2:end-2)) + f(end));


%%
clear; close; clc;

rng(1);         
x = randn(10000,1);

mean(x)
std(x)

histogram(x,40,'Normalization','pdf');
hold on
fplot(@(z) normpdf(z,0,1),[-4,4],'LineWidth',2)

%%
clear; close; clc;

R = 2000; n = 100;
beta1 = zeros(R,1);

for r = 1:R
    x = randn(n,1);
    u = randn(n,1);
    y = 1 + 0.5*x + u;
    bhat = regress(y,[ones(n,1) x]);
    beta1(r) = bhat(2);
end

histogram(beta1,40,'Normalization','pdf'); hold on
mu = mean(beta1); 
se = std(beta1);
fplot(@(b) normpdf(b,0.5,se),[0,1],'LineWidth',2)