clear all;

global delta alpha A ksi 
% beta = 0.985;
% delta = 0.025;
% alpha = 0.35;
% A = 1;
% omega_bar = 0.831;
% chi = 6.43;
% ksi = 0.5;
% eta = 0.164; 
% sigma = 1/eta;
% eps_min = 1-eta;
% ppi = 0.00324;
% B = 0.618;
beta = 0.9848;
delta = 0.025;
alpha = 0.35;
A = 1;
omega_bar = 0.76;
chi = 6.43;
ksi = 0.5;
eta = 0.1639; 
sigma = 1/eta;
eps_min = 1-eta;
ppi = 0.00305;
gamma = (1-delta)/(1-ppi);
B = 0.618;
pam0 = [beta chi omega_bar ppi eta B];

[pam,fval,exitflag,output]=fminunc(@(pam)solvepam(pam,delta,alpha,A,ksi),pam0);

% options = optimset('TolFun',1e-5,'MaxFunEvals',1e4);
% [pam,fval,exitflag,output]=fsolve(@(pam)solvepam(pam,delta,alpha,A,ksi),pam0,options);
solvepam(pam,delta,alpha,A,ksi)

function distance = solvepam(pam,delta,alpha,A,ksi)
    beta = pam(1);
    chi = pam(2);
    omega_bar = pam(3);
    ppi = pam(4);
    eta = pam(5);
    B = pam(6);
    sigma=1/eta;
    eps_min = 1-eta;
    NN = 500; % number of grids

    qk_ss = 0.6677;

    qk_lb = 0.93*qk_ss;
    qk_ub = 1.1*qk_ss;
    
    qk = linspace(qk_lb,qk_ub, NN);
    
    eps2 = 1./qk;
    
    theta = ((eps2./eps_min).^(-sigma)).*(1-delta)./(ppi+((eps2./eps_min).^(-sigma)).*(1-delta));
    
    Phi_K = qk./theta;
    
    eps1 = 1./Phi_K;
    
    zeta = (eps1./eps_min).^(-sigma);
    
    qb = beta.*(1+zeta./(sigma-1));
    
    R = 1./qb;
    
    rk = (1-beta.*(1-delta).*((1+(theta.^sigma).*zeta./(sigma-1)).*omega_bar+(1-omega_bar))).*Phi_K.*R-omega_bar.*ppi.*qk;
    
    B_K = delta./(zeta.^(1-1./sigma))-(rk+omega_bar.*ppi.*qk)-(1-delta).*omega_bar.*qk.*(theta.^(sigma-1));
    
    I_K = (B_K+rk+omega_bar.*ppi.*qk).*zeta+(1-delta).*omega_bar.*qk.*(theta.^sigma).*zeta;
    
    Y_K = rk./alpha;
    
    C_K = Y_K-I_K-(1-qb).*B_K;
    
    N = (((1-alpha)./chi).*(Y_K./C_K)).^(1./(1+ksi));
    
    K = ((alpha.*A./rk).^(1./(1-alpha))).*N;
    
    omega = (ppi./(ppi+1-delta)+(1-delta).*(theta.^sigma)./(ppi+1-delta)).*omega_bar;
    
    Y = Y_K.*K;
    
    I = I_K.*K;
    
    C = C_K.*K;
    
    BB = B_K.*K;
    
    % welfare
    V = (log(C)-chi.*(N.^(1+ksi))./(1+ksi))./(1-beta);
    
    B_Y = B_K./Y_K;
    
    I_Y = I_K./Y_K;
    
    w = (1-alpha).*Y./N;
    
    Rk = R./theta;
    
    S = (ppi+(1-delta).*(theta.^sigma).*zeta).*K;
    
    spread = Rk-R;
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % TO SOLVE AND SAVE THE STEADY STATE  %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    [error,a] = min(abs(BB-B*ones(1,NN)));
    
    qk_ss = qk(a);
    Phi_k_ss = Phi_K(a);
    qb_ss = qb(a);
    B_ss = BB(a);
    C_ss = C(a);
    I_ss = I(a);
    K_ss = K(a);
    Y_ss = Y(a);
    N_ss = N(a);
    R_ss = R(a);
    rk_ss = rk(a);
    zeta_ss = zeta(a);
    omega_ss = omega(a);
    eps1_ss = eps1(a);
    eps2_ss = eps2(a);
    theta_ss = theta(a);
    BY_ss = B_Y(a);
    IY_ss = I_Y(a);
    w_ss = w(a);
    Lambda_ss = 1/C_ss;
    Rk_ss = Rk(a);
    spread_ss = spread(a);
%     l_ppi_ss = R_ss*beta*(1-theta_ss^(1/eta))*(eta/(1-eta))*((eps1_ss/eps_min)^(-1/eta))*(1-delta)*omega_bar-R_ss*beta*(1+(eta/(1-eta))*((eps1_ss/eps_min)^(-1/eta)))*omega_bar*ppi*qk_ss/Phi_k_ss;
%     l_omega_ss = R_ss*beta*((eta/(1-eta))*((eps1_ss/eps_min)^(-1/eta)))*(1-delta)*(1-omega_bar);
    l_ppi_ss = R_ss*beta*(1-theta_ss^(1/eta))*(eta/(1-eta))*((eps1_ss/eps_min)^(-1/eta))*(1-delta)-R_ss*beta*(1+(eta/(1-eta))*((eps1_ss/eps_min)^(-1/eta)))*ppi*qk_ss/Phi_k_ss;
    l_omega_ss = R_ss*beta*(1+(eta/(1-eta))*((eps1_ss/eps_min)^(-1/eta)))*ppi*(1-omega_bar)*qk_ss/Phi_k_ss+R_ss*beta*((eta/(1-eta))*((eps1_ss/eps_min)^(-1/eta)))*(1-delta)*(1-omega_bar)-R_ss*beta*(1-theta_ss^(1/eta))*(eta/(1-eta))*((eps1_ss/eps_min)^(-1/eta))*(1-omega_bar)*(1-delta);

    dis1 = BY_ss-0.6;
    dis2 = IY_ss-0.2;
%     dis3 = 100*(R_ss-1.0055);
    dis3 = R_ss-1.0055;
    dis4 = N_ss-0.25;
    dis5 = 100*l_ppi_ss - 0.2595;
    dis6 = 100*l_omega_ss - 0.1608;

    %distance = [dis1,dis2,dis3,dis4,dis5,dis6];
    distance = sqrt(dis1^2+dis2^2+dis3^2+dis4^2+dis5^2+dis6^2);

end