var Y K H r w cf cd tot h R q i xd phi PHI f g n ke ce B cc Cc Ce C I b wr A omegab rif TB xf G;

varexo eA eP;

parameters eps theta psi alpha zeta beta delta  kapd kapf gamma eta S mu M rhoA rhoP StdeA StdeP wr_st sq kap;      


         alpha    = 0.36;  
		 zeta     = 1-alpha; 
         eps = 1;
         theta = 0.5;
         psi  = 2;
		 beta     = 0.99; 
		 delta    = 0.02; 
         kapd     =2;
         kapf     =2;
        gamma     =0.9474;
        eta       =0.1;
          S =  0.207;
          M        = -.5*S^2; 
          mu  = 0.25;
          wr_st = 0.01;
        StdeA = 0.01; 
         rhoA = 0.99;
         rhoP= 0.8;
       StdeP = 0.01;
          sq = 0.004;
          kap = 0.1;
    



model;
//production
    G = A*(K(-1)^alpha)*(H^zeta);
    r = alpha*A*(K(-1)^(alpha-1))*(H^zeta);// MPK
    w = zeta*A*(K(-1)^alpha)*(H^(zeta-1)); //MPH

//consumer FOC
(theta/(1-theta))*(cf/cd)=1/tot;
(psi/((1-h)*theta))*(((cd^theta)*(cf^(1-theta)))^eps)*(cf/cd)^(theta-1)=w;
(R(+1)/beta) = ((((cd^theta)*(cf^(1-theta)))/((cd(+1)^theta)*(cf(+1)^(1-theta))))^eps)*((cf(+1)/cf)*(cd/cd(+1)))^(1-theta);
(1/q)*((1-delta)*q(+1)+r(+1))= 1/R(+1);

//investment
i = kapd*xd;
xf = (kapd/kapf)*xd;

//optimal contract
    q*kapd*(1 - mu*PHI - (mu*phi*f)/(1-PHI)) = 1+(kapd/kapf)*tot; //equation7and11
    xd = (1/(1+(kapd/kapf)*tot-q*kapd*g)) * n; //equation8 and 12

//entrepreneur
    n =(ke(-1))*(q*(1-delta)+r); //equation9 and 15
    (q/tot) =  (beta*gamma) * ((r(+1) + q(+1)*(1-delta))*rif/tot(+1))  * ( (q(+1)*f(+1))/(1-q(+1)*g(+1)) ); //equation 11
    ke = i*f - (ce*tot)/q; //equation 10

   rif = q*f*i/n;

//Capital accumulation
    K = (1-delta)*K(-1) + I*(1-mu*PHI);//equation5



//market clearing
     r*K(-1)+w*h+q*(1-delta)*K(-1)+(b*R)-b(-1)-cd-tot*cf-q*K=0; //equation19 check
    (1-eta)*(cd+tot*cf) + eta*ce + eta*i+(1-eta)*b(-1) = Y; //MC4
    H = (1-eta)*h; //MC1
    cc = cd+tot*cf;
    Cc = (1-eta)*(cd+tot*cf);       
    Ce = eta*ce;            
    C = (1-eta)*cc + eta*ce; 
    I = eta*i;  
     B(-1) = (1-eta)*b(-1);
    TB = G -(1-eta)*(tot*cf+cd)-eta*(xd+tot*xf);

//lognormal                 
    PHI = normcdf((log(omegab)-M)/S);  
    phi = normpdf((log(omegab)-M)/S) / (omegab*S);                                    
    g   = normcdf((log(omegab)-M)/S - S) - PHI*mu + (1-PHI)*omegab; 
    f   = 1 - mu*PHI - g;    

// domestic interest rate equation
R = 1/(1+wr);                               // R is discount facot and wr is real interest rate
R(+1)=(1/(1+wr_st))*exp(-sq*B(-1)-kap*(B-B(-1)));   // domestic real interest rate pinned down by the exogensour world interest rate wr_st and the surplus or deficit level

//shocks
    log(A)    = rhoA*log(A(-1)) + StdeA*eA;
    log(tot) =  rhoP*log(tot(-1))+StdeP*eP;
   
end;


initval;
    omegab  = 0.603892;
    PHI     = normcdf((log(omegab)-M)/S);
    phi     = normpdf((log(omegab)-M)/S) / (omegab*S);
    g       = normcdf((log(omegab)-M)/S - S) - PHI*mu + (1-PHI)*omegab;
    f       = 1-mu*PHI- g;
    q       = 1/(1-mu*0.01+(gamma-1)*f);
    r       = q*((1-beta*(1-delta))/beta);
    H       = .3;
    h      = .3/(1-eta);
    K       = (alpha/r)^(1/(1-alpha))*(H^(zeta/(1-alpha)));
    Y       = (K^alpha)*(H^zeta);
    i       = (delta/(eta*(1-mu*0.01)))*K;
    n       = (1-g*q)*i;
    ke      = (beta/q)*(eta*n-(1-alpha-zeta)*Y);
    ce      = q*(f*i-(ke/eta));
    cc      = (Y - eta*ce - eta*i)/(1-eta);
    cf      = 0.45;
    cd      =0.45;
    A       = 1;
    tot     = 1;
      R     = 0.99;
      wr    = (1/R) -1;
end;

steady;
check;

shocks;
    var eA = 1;
    var eP = 1;
end;

stoch_simul(order=1,irf=24);