function parameter=calibration
clear;

R=1+3.87/400;   
PI=1+2/400;     
BETA=1/R*PI;    
DELTA=0.025;
THETA_S=0.7;      
THETA_P=0.4;
ALPHA=0.5;
As=1;
Ap=1.42;
ETA=1;
GAMMA=0.1;
TAU_S=0.5;
TAU_P=0.4;
SIGMA=3; 
qk=1;
F_xi=1.83/400;  
MU=0.15;
ZETA=0;
OMEGA_I=2.5; 
RL=1+2/400;
Rp=1+4.8/400;
Rs=1+4.2/400;
omegas=1-0.7239; 
[lnr,fval,nflagr]=fsolve(@(x) solve_calim(x(1),omegas,F_xi),1/0.4,optimset('TolFun',1e-10));
KM=lnr(1);OMEGAM=1-1/KM;
int_xi=KM/(KM-1)*(OMEGAM-OMEGAM^KM*omegas^(1-KM));
%F_xi=1-(OMEGAM/omegas)^KM;
Xb=1-F_xi+(1-MU)*int_xi/omegas;
PHI=0.5;
BETA_B=0.9952;
lambda_c=1-BETA_B*R/PI;

rb=(-lambda_c*(1-GAMMA*ZETA)+1)/BETA_B*PI/Xb;

BETA_SB=1/(1-F_xi)/rb*PI;
Rp_bar=Rs; 
EPSILON=10;
pw=(EPSILON-1)/EPSILON;
BETA_E=BETA;
rk_s=qk*(1-THETA_S/Rs*PI)/BETA_E-(1-DELTA)+THETA_S*qk;
XI=0.21;
rk_p=qk*(1-THETA_P/Rp*PI-XI)/BETA_E-(1-DELTA)+THETA_P*qk+XI*Rp_bar/PI;
[lir,fval,cflagr]=fsolve(@(x) solve_cali2(x(1),x(2),SIGMA,PHI,ALPHA,rk_s,rk_p,As,Ap,pw),[0.5,0.5],optimset('TolFun',1e-10));
ps=lir(1);pp=lir(2);
ys_y=(ps/pw)^(-SIGMA)*PHI^(SIGMA);
yp_y=(pp/pw)^(-SIGMA)*(1-PHI)^(SIGMA);
ks_hs=(rk_s/ps/ALPHA/As)^(1/(ALPHA-1));
kp_hp=(rk_p/pp/ALPHA/Ap)^(1/(ALPHA-1));

[ltr,fval,tflagr]=fsolve(@(x) solve_cali3(x(1),ks_hs,kp_hp,As,Ap,ALPHA,ys_y,yp_y),1/4,optimset('TolFun',1e-10));
hs=ltr(1);
hp=1/3-hs;
ks=ks_hs*hs;
kp=kp_hp*hp;
ys=As*ks^ALPHA*hs^(1-ALPHA);
yp=Ap*kp^ALPHA*hp^(1-ALPHA);
y=ys/ys_y;
wage=(1-ALPHA)*ps*ys/hs;
ls=THETA_S*qk*ks/Rs*PI;
lp=THETA_P*qk*kp/Rp*PI;
lp_bar=XI*kp;
b=lp_bar*Rp_bar*omegas/rb;
ALPHA_S=(BETA_B*(Rs-RL*TAU_S)/PI-(1-TAU_S)+lambda_c*(1-GAMMA))/ls;
ALPHA_P=(BETA_B*(Rp-RL*TAU_P)/PI-(1-TAU_P)+lambda_c*(1-GAMMA))/lp;
ALPHA_SB=(BETA_SB*Rp_bar/PI*(1-int_xi)-1)/lp_bar;
DWL=int_xi*MU*lp_bar*Rp_bar;
d=lp+ls+b-GAMMA*(ls+lp+ZETA*b);
PI_SB=((1-int_xi)*lp_bar*Rp_bar-(1-F_xi)*rb*b)/PI-lp_bar+b-ALPHA_SB/2*lp_bar^2;
PI_CB=(Rp*lp+Rs*ls+rb*b*Xb-R*d-RL*(TAU_P*lp+TAU_S*ls))/PI-ALPHA_P/2*lp^2-ALPHA_S/2*ls^2-(1-TAU_S)*ls-(1-TAU_P)*lp-b+d;

DIV_P=(-DELTA)*kp+pp*yp-wage*hp+(1-Rp/PI)*lp+(1-Rp_bar/PI)*lp_bar;
DIV_S=(-DELTA)*ks+ps*ys-wage*hs+(1-Rs/PI)*ls;
k=kp+ks;
i=DELTA*k;
ch=y-DWL-i;
KAPPA=wage*3^ETA/ch;
m=d/95*5;
bg=((TAU_P*lp+TAU_S*ls)*(RL/PI-1)+(1-1/PI)*m)/(1-R/PI);
lambda_h=1/ch;
KAPPA_M=m*(lambda_h-lambda_h/PI*BETA);
OMEGA_P=22;
ZETAL=1;

parameter=[BETA DELTA THETA_S THETA_P ALPHA As Ap ETA GAMMA TAU_S TAU_P SIGMA OMEGA_I MU ZETA EPSILON PI RL KM BETA_B BETA_E PHI ALPHA_S ALPHA_P OMEGA_P KAPPA XI KAPPA_M BETA_SB ALPHA_SB ZETAL OMEGAM];