% This function takes as input the wage parameters and generates:
% (1) A wage grid (F_t-1) by age, that is growing by the rate of sigma(v)
% for each earner
% (2) Grids for v1 v2 u1 u2 which are constant by age. Note that the
% correlation structure between husband and wife implies these are grids
% with probability on each combination.

function wage_grids_outputs = wage_grids(parameters)
TT      = parameters.TT     ;
div     = parameters.div    ;
minage  = parameters.minage ;
wpoints = parameters.wpoints;
vpoints = parameters.vpoints;
upoints = parameters.upoints;
theta   = parameters.theta  ;
maxw    = parameters.maxw;

%% Set Parameters 
sd_u1   = sqrt(theta(1));
sd_u2   = sqrt(theta(2));
sd_v1   = sqrt(theta(3));
sd_v2   = sqrt(theta(4));
cov_u1u2= theta(5);
cov_v1v2= theta(6);
sd_w0_1 = sqrt(theta(7));
sd_w0_2 = sqrt(theta(8));
detw_ag_1 = theta(9);
detw_ag_2 = theta(10);
detw_ag2_1= theta(11);
detw_ag2_2= theta(12);
detw_c_1  = theta(13);
detw_c_2  = theta(14);

%% Generate the deterministic path 
% Generate the annual deterministic path
TTT         = TT/div;
maxage      = minage + TTT;
agegrid     = linspace(minage,maxage,maxage-minage+1);
age_mat     = [ones(size(agegrid))',agegrid',(agegrid.^2)'];
beta        = [detw_c_1; detw_ag_1; detw_ag2_1];
log_w_det_an_1 = age_mat*beta; 
beta        = [detw_c_2; detw_ag_2; detw_ag2_2];
log_w_det_an_2 = age_mat*beta; 

% Apply the subdivision
log_w_det_1 = zeros(TT+1,1);
for tau=1:TTT
    startw = log_w_det_an_1(tau);
    step   = (log_w_det_an_1(tau+1)-log_w_det_an_1(tau))/div;
    endw   = log_w_det_an_1(tau+1)-step;
    log_w_det_1((tau-1)*div+1:tau*div) = linspace(startw,endw,div);
end
log_w_det_1(TT+1)=log_w_det_an_1(TTT+1);

log_w_det_2 = zeros(TT+1,1);
for tau=1:TTT
    startw = log_w_det_an_2(tau);
    step   = (log_w_det_an_2(tau+1)-log_w_det_an_2(tau))/div;
    endw   = log_w_det_an_2(tau+1)-step;
    log_w_det_2((tau-1)*div+1:tau*div) = linspace(startw,endw,div);
end
log_w_det_2(TT+1)=log_w_det_an_2(TTT+1);

%% Factorize the shocks 
sig_u = [sd_u1^2 cov_u1u2;cov_u1u2 sd_u2^2];
sig_v = [sd_v1^2 cov_v1v2;cov_v1v2 sd_v2^2];
[Lu,Du] = ldl(sig_u);
[Lv,Dv] = ldl(sig_v);

gamma_u   = Lu(2,1);
gamma_v   = Lv(2,1);
sig_u2_hat = Du(2,2);
sig_v2_hat = Dv(2,2);

%% Generate grid for the permanent and transitory shocks
[v1_xgh, v1_wgh] = hermquad(vpoints);
v1_xgh = flip(sqrt(2)*sd_v1*v1_xgh);
v1_wgh = flip((1/sqrt(pi))*v1_wgh);

[v2_xgh, v2_wgh] = hermquad(vpoints);
v2_xgh = flip(sqrt(2)*sqrt(sig_v2_hat)*v2_xgh);
v2_wgh = flip((1/sqrt(pi))*v2_wgh);

[u1_xgh, u1_wgh] = hermquad(upoints);
u1_xgh = flip(sqrt(2)*sd_u1*u1_xgh);
u1_wgh = flip((1/sqrt(pi))*u1_wgh);

[u2_xgh, u2_wgh] = hermquad(upoints);
u2_xgh = flip(sqrt(2)*sqrt(sig_u2_hat)*u2_xgh);
u2_wgh = flip((1/sqrt(pi))*u2_wgh);

%% Grid by age 
v1_grid = [v1_xgh*ones(1,TT)];
v2_grid = [v2_xgh*ones(1,TT)];
u1_grid = [u1_xgh*ones(1,TT)];
u2_grid = [u2_xgh*ones(1,TT)];

% adjust for idle state assuming shocks hit only once a year
zers    = zeros(floor(vpoints/2),TT);
probdiv = ((div-1)/div)*ones(1,TT);
v1_prob = [(1/div)*v1_wgh*ones(1,TT)] + [zers; probdiv; zers];
v2_prob = [(1/div)*v2_wgh*ones(1,TT)] + [zers; probdiv; zers];

zers    = zeros(floor(upoints/2),TT);
u1_prob = [(1/div)*u1_wgh*ones(1,TT)] + [zers; probdiv; zers];
u2_prob = [(1/div)*u2_wgh*ones(1,TT)] + [zers; probdiv; zers];

%% Generate equally spaced grid for the log wage by age (mu is log_w_det) SD growing according to the permanent component
% grid is generated as normal 0,1 with maxw S.Ds for each direction
xgh = linspace(-maxw,maxw,wpoints);

var_age_1    = zeros(size(log_w_det_1));
var_age_2    = var_age_1;
var_age_1(1) = (sd_w0_1)^2 - (sd_u1)^2;
var_age_2(1) = (sd_w0_2)^2 - (sd_u2)^2;

for i=2:TT
   var_age_1(i) = var_age_1(i-1) + (sd_v1^2)/div; 
   var_age_2(i) = var_age_2(i-1) + (sd_v2^2)/div; 
end
sd_age_1 = sqrt(var_age_1);
sd_age_2 = sqrt(var_age_2);

lw_grid_1 = zeros(wpoints,TT);
lw_grid_2 = zeros(wpoints,TT);
for i=1:TT 
   lw_grid_1(:,i) = sd_age_1(i)*xgh + log_w_det_1(i);
   lw_grid_2(:,i) = sd_age_2(i)*xgh + log_w_det_2(i);
end

%% Store results
wage_grids_outputs.v1_grid = v1_grid;
wage_grids_outputs.v2_grid = v2_grid;
wage_grids_outputs.u1_grid = u1_grid;
wage_grids_outputs.u2_grid = u2_grid;
wage_grids_outputs.v1_prob = v1_prob;
wage_grids_outputs.v2_prob = v2_prob;
wage_grids_outputs.u1_prob = u1_prob;
wage_grids_outputs.u2_prob = u2_prob;
wage_grids_outputs.lw_grid_1 = lw_grid_1;
wage_grids_outputs.lw_grid_2 = lw_grid_2;
wage_grids_outputs.gamma_u    = gamma_u; 
wage_grids_outputs.gamma_v    = gamma_v;   
wage_grids_outputs.sig_u2_hat = sig_u2_hat;
wage_grids_outputs.sig_v2_hat = sig_v2_hat;
wage_grids_outputs.log_w_det_1=log_w_det_1;
wage_grids_outputs.log_w_det_2=log_w_det_2;


