%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Learning model solver: solve the learning model
%	  -  input: the PLM beliefs and innovations and lagged R, y, k , gamma
%	  -  output: endogenous variables under learning
%	  -  enter the matrix and invert it to solve 
%
%      xw : trend belief in w equation
%      xR : trend belief in R equation
%      gam: current shock
%      lagk : lagged k  (means the state)
% ini_con: shockmat(i,t) khat1 ahat1 bhat1


%%%
% The problem is JJ(2,1) is non-zero
%


function [y_vec1, fore_dis]  = Learningsol_RE_gshock(xw, xR, xk, xb, ini_con)
%function [y_vec1, coeff_mat] = Learningsol(xw, xR, xk, xb, ini_con)
global delta alpha gamma eps_H eta sigma tauk tauH tauC R_bar beta_tilda gamma1 gamma2 lam1 lam2 ...
    theta ik_ratio  cy_ratio iy_ratio gy_ratio ky_ratio ay_ratio by_ratio gt at bt rt ...
   R_tilda psi rk wage output  hours cons invst ...
    caput cap gamma_x coeff_vec_k_RE coeff_vec_s_RE  epsZ F  epsc epsw kai ini_wR
% 


gam = ini_con(1); 
khat1 = ini_con(2); 
%ahat1 = ini_con(3); 
bhat1 = ini_con(3); 
gshock = ini_con(4);
% display('gov shock')
% gshock

w0vec = [xk(1); xb(1)];
Amatrix_AL = [xk(2) xk(3); xb(2) xb(3)]; 
Cmatrix_AL = [xw(2) xw(3) ; xR(2)  xR(3)];
w0tilda = [xw(1) ; xR(1)]; 

ww = [(- khat1 + log(gamma) +gam); (- bhat1 + log(gamma) +gam)]; 

IA = inv(eye(2)-Amatrix_AL); 
IBA= inv(eye(2)-beta_tilda*Amatrix_AL); 

% display('(eye(2)-Amatrix_AL)\w0vec==log(gamma)')
% 
% (eye(2)-Amatrix_AL)\w0vec == [log(gamma);log(gamma)]

TT = beta_tilda/(1-beta_tilda)*Cmatrix_AL*IBA;

% GG= beta_tilda/(1-beta_tilda)^2*(w0tilda + Cmatrix_AL*inv(eye(2)-Amatrix_AL)*w0vec) ...
%     +TT * IA*((eye(2)-Amatrix_AL)*ww-w0vec);

GG= beta_tilda/(1-beta_tilda)^2*(w0tilda + Cmatrix_AL*inv(eye(2)-Amatrix_AL)*w0vec) ...
    +TT * (ww-IA*w0vec);

% GG= beta_tilda/(1-beta_tilda)^2*(w0tilda + Cmatrix_AL*IA*w0vec) ...
%     +TT * IA*((eye(2)-Amatrix_AL)*ww-w0vec);

% display('(eye(2)-Amatrix_AL)*ww-w0vec')
% (eye(2)-Amatrix_AL)*ww-w0vec
% display('ww')
% ww
% display('IA')
% IA
% display('w0vec')
% w0vec

% display('ww(1)==ww(2)')
% ww(1)==ww(2)
% display('Amatrix_AL')
% Amatrix_AL
% display('(eye(2)-Amatrix_AL)*ww-w0vec')
% ttt= (eye(2)-Amatrix_AL)*ww-w0vec;
% ttt(1)
% ttt(2)
% 
% display('-Amatrix_AL(2,1)*ww(1)+(1- Amatrix_AL(2,2))*ww(2) -w0vec(2)    ')
% -ww(1)*Amatrix_AL(2,1)+ww(2)*(1- Amatrix_AL(2,2)) -w0vec(2)

% display('TT * IA*((eye(2)-Amatrix_AL)*ww-w0vec)')
% 
% TT * IA*((eye(2)-Amatrix_AL)*ww-w0vec)


% display('beta_tilda/(1-beta_tilda)^2*(w0tilda + Cmatrix_AL*IA*w0vec)')
% 
% beta_tilda/(1-beta_tilda)^2*(w0tilda + Cmatrix_AL*IA*w0vec)


% Del0 = -(1-kai)*(1-beta_tilda)/epsc*beta_tilda^(-1)*(by_ratio/ay_ratio*bhat1+ky_ratio/ay_ratio*khat1-gam)...
%     -gamma1*GG(2)-gamma2*(GG(1)-beta_tilda/(1-beta_tilda)^2*log(gamma)); 

Del0 = (1-kai)*(1-beta_tilda)/epsc*beta_tilda^(-1)*(by_ratio/ay_ratio*bhat1+ky_ratio/ay_ratio*khat1-gam)...
    +gamma1*(GG(2)-beta_tilda/(1-beta_tilda)^2*ini_wR(2))+gamma2*(GG(1)-beta_tilda/(1-beta_tilda)^2*ini_wR(1)); 

%Del0 =0 ; 

% display('TT')
% 
% TT
% 
% display('GG')
% GG
% display(' GG(2)-beta_tilda/(1-beta_tilda)^2*ini_wR(2)')
% GG(2)-beta_tilda/(1-beta_tilda)^2*ini_wR(2)
% display('-(1-kai)*(1-beta_tilda)/epsc*beta_tilda^(-1)*(by_ratio/ay_ratio*bhat1+ky_ratio/ay_ratio*khat1-gam)')
% -(1-kai)*(1-beta_tilda)/epsc*beta_tilda^(-1)*(by_ratio/ay_ratio*bhat1+ky_ratio/ay_ratio*khat1-gam)
% display('gamma1*(GG(2)-beta_tilda/(1-beta_tilda)^2*ini_wR(2))')
% 
% gamma1*(GG(2)-beta_tilda/(1-beta_tilda)^2*ini_wR(2))
% 
% display('gamma2*(GG(1)-ini_wR(1))')
% gamma2*(GG(1)-beta_tilda/(1-beta_tilda)^2*ini_wR(1))
% 
% display('beta_tilda/(1-beta_tilda)^2*ini_wR(1)')
% beta_tilda/(1-beta_tilda)^2*ini_wR(1)
% display('bhat1')
% bhat1
% display('khat1')
% khat1
% 
% display('gam')
% gam
% 
% display('Del0')
% Del0
% display('GG(1)')
% GG(1)
% display('GG(1)-beta_tilda/(1-beta_tilda)^2*ini_wR(1)')
% GG(1)-beta_tilda/(1-beta_tilda)^2*ini_wR(1)
% 
% display('GG(2)-ini_wR(2)')
% GG(2)-ini_wR(2)

Del1 = sigma^(-1) * psi*(1-sigma);
Del2 = 1;
Del3 = -gamma2/(1-beta_tilda); 
Del4 = -(1-kai)*(1-beta_tilda)/epsc*R_tilda -gamma1*beta_tilda/(1-beta_tilda); 
Del5 = -gamma1*TT(2,1) -gamma2*TT(1,1);
Del6 = -gamma1*TT(2,2) -gamma2*TT(1,2);


FF = zeros(11,11);

FF(1,wage) = 1;
FF(1,hours) =-(eps_H-psi*(sigma-1)/sigma);
FF(1,cons) = -1; 

FF(2, wage) = -1; 
FF(2, output) = 1;
FF(2, hours) = -1; 

FF(3, rk) = -1;
FF(3, output) = 1; 
FF(3, caput) = -1;

FF(4,output) = -1; 
FF(4,cons) = cy_ratio; 
FF(4, invst) = iy_ratio;
FF(4,gt) = gy_ratio; 

FF(5, output) = -1; 
FF(5, hours) = 1-alpha; 
FF(5, caput) = alpha; 

FF(6, invst) = ik_ratio ;
FF(6, caput) = -delta*theta/gamma; 
FF(6, cap) = -1; 

FF(7, rk) = 1/(1-theta);
FF(7, caput) = 1; 

FF(8, rk) = Del4; 
FF(8, wage) = Del3; 
FF(8, hours) = Del1; 
FF(8, cons) = Del2;
FF(8, cap) = Del5;
FF(8, bt) = Del6;

FF(9, output) = gy_ratio+lam1 - lam2 *by_ratio/gamma;
FF(9, gt) = -gy_ratio; 

FF(10, rk) = - tauk*ky_ratio*theta*delta/gamma + by_ratio*R_bar/gamma*theta*delta*(1-tauk)/R_bar; 
FF(10, wage) = -(1-alpha)*tauH;
FF(10, hours) = -(1-alpha)*tauH ;
FF(10, gt) = gy_ratio;
FF(10, bt) = -by_ratio; 
FF(10,cons) = -cy_ratio*tauC; 


FF(11, gamma_x) = 1; 

% display('FF')
% FF
% display('max(abs(eig(FF)))')
% max(abs(eig(FF)))
% 
% display('conditional number')
% 
% cond(FF)

G = zeros(11,1);

G(3) = khat1-gam;

G(5) = -alpha*G(output); 

G(6) = -(1-delta)/gamma*G(output); 

G(8) = Del0; 
G(9) = -lam2  * by_ratio / gamma * (bhat1-gam) -gy_ratio*gshock;  
G(10) =  -by_ratio *R_bar / gamma *bhat1  ...
    -( (R_bar-1)/(1-tauk) *tauk *ky_ratio/gamma -by_ratio*R_bar/gamma)*gam...
    + (R_bar-1)/(1-tauk)*tauk/gamma*ky_ratio * khat1 ;

G(11) = gam; 

% display('G value')
% G

y_vec1 = FF\G; 

wage1 = y_vec1(2);
cap1 = y_vec1(9); 
bt1 = y_vec1(11); 
rk1 = y_vec1(1) ; 

swe = beta_tilda/ (1-beta_tilda) * wage1 +GG(1)+TT(1,1)*cap1...
    +TT(1,2)*bt1-beta_tilda/(1-beta_tilda)^2*log(gamma); 
swr = beta_tilda/(1-beta_tilda) * rk1 + GG(2) +TT(2,1)* cap1 + TT(2,2) * bt1; 

% display('GG(1)-...')
% 
% GG(1) - beta_tilda/(1-beta_tilda)^2*log(gamma)

fore_dis =  [swe; swr]; 
















































