function [PP,QQ,RR,SS,NN,PROBLEM]=modelehl(sigma,thetap,thetaw,gam,phi,alpha,rhor,gy,gpi,rhoa,rhog);

% Fixed Parameters

beta = 0.99;
epsilon=6;
delta=0.36;


% Coefficients

kappap=(1-thetap)*(1-thetap*beta)*(1-delta)/(thetap*(1+delta*(epsilon-1)));
kappaw=(1-thetaw)*(1-thetaw*beta)/(thetaw*(1+phi*gam));


% 3) 0 = -mrs(t) + c(t) / sigma + gam * n(t) - g(t)
% 4) 0 = -i(t) + rhoi * i(t-1) + gammapi * (1-rhoi) * dp(t) + gammay * (1-rhoi) * y(t) + ms(t)
% 5) 0 = -y(t) + c(t)
% 6) 0 = mrs(t) - x(t)
% 7) 0 = E_t [ -sigma * i(t) + sigma * dp(t+1) - sigma * g(t+1) + sigma * g(t) - c(t) + c(t+1)]
% 8) 0 = E_t [ kappap * mc(t) + kappap * mu(t) - dp(t) + beta * dp(t+1)]
% 9) a(t+1) = rhoa * a(t) + epsilona(t+1)
% 10) ms(t+1) = epsilonm(t+1)
% 11) mu(t+1) = epsilonp(t+1)
% 12) g(t+1) = rhog * g(t) + epsilong(t+1)


% The equations are, conveniently ordered:

% 1) 0 = -y(t) + a(t) + (1-delta) * n(t) 
% 2) 0 = y(t) - n(t) + mc(t) - x(t)
% 3) 0 = -mrs(t) + c(t) / sigma + gam * n(t) - g(t)
% 4) 0 = -i(t) + rhoi * i(t-1) + gammapi * (1-rhoi) * dp(t) + gammay * (1-rhoi) * y(t) + ms(t)
% 5) 0 = -y(t) + c(t)
% 6) 0 = -x(t) + x(t-1) + dw(t) - dp(t)
% 7) 0 = E_t [ -sigma * i(t) + sigma * dp(t+1) - sigma * g(t+1) + sigma * g(t) - c(t) + c(t+1)]
% 8) 0 = E_t [ kappap * mc(t) + kappap * mu(t) - dp(t) + beta * dp(t+1)]
% 9) 0 = E_t [ kappaw * mrs(t) - kappaw * x(t) - dw(t) + beta * dw(t+1)]
% 10) a(t+1) = rhoa * a(t) + epsilona(t+1)
% 11) ms(t+1) = epsilonm(t+1)
% 12) mu(t+1) = epsilonp(t+1)
% 13) g(t+1) = rhog * g(t) + epsilong(t+1)
%
% Endogenous state variables "x(t)": x(t), i(t), dp(t), dw(t), y(t)
% Endogenous other variables "y(t)": n(t), mc(t), mrs(t), c(t)
% Exogenous state variables  "z(t)": a(t), ms(t), mu(t), g(t).
% Switch to that notation.  Find matrices for format
% 0 = AA x(t) + BB x(t-1) + CC y(t) + DD z(t)
% 0 = E_t [ FF x(t+1) + GG x(t) + HH x(t-1) + JJ y(t+1) + KK y(t) + LL z(t+1) + MM z(t)]
% z(t+1) = NN z(t) + epsilon(t+1) with E_t [ epsilon(t+1) ] = 0,

%Order:x,i,dp,dw,y

AA = [ 0,  0,  0,  0, -1,
      -1,  0,  0,  0,  1,
       0,  0,  0,  0,  0,
       0,  -1, gpi*(1-rhor), 0, gy*(1-rhor),
      -1,  0, -1,  1,  0,
       0,  0,  0,  0, -1];
       
BB = [ 0,  0,  0,  0, 0,
       0,  0,  0,  0, 0,
       0,  0,  0,  0, 0,
       0, rhor, 0, 0, 0,
       1,  0,  0,  0, 0,
       0,  0,  0,  0,  0];
       
%Order:  n	     mc	  mrs	  c      
CC = [  1-delta,  0,    0,   0,    		% Equ. 1)
		  -1,       1,    0,   0,    		% Equ. 2)
	    gam,       0,   -1,   1/sigma   % Equ. 3)
         0,       0,    0,   0,        % Equ. 4)
         0,       0,    0,   0,        % Equ. 5)
         0,       0,    0,   1];      % Equ. 6)

      
DD = [ 1, 0, 0, 0,
       0, 0, 0, 0, 
       0, 0, 0, -1,
       0, 1, 0, 0, 
       0, 0, 0, 0,
       0, 0, 0, 0];
    
FF = [ 0, 0, sigma, 0, 0,
       0, 0, beta, 0, 0,
       0, 0, 0, beta, 0];
    
GG = [ 0, -sigma, 0, 0, 0,
       0, 0, -1, 0, 0,
      -kappaw, 0, 0, -1, 0];
    
HH = [ 0, 0, 0, 0, 0,
       0, 0, 0, 0, 0,
       0, 0, 0, 0, 0 ];
    
JJ = [ 0,  0,  0,  1,  
		 0,  0,  0,  0, 
		 0,  0,  0,  0];

KK = [  0,  0,  0,  -1,  
		  0, kappap,  0,  0,    
		  0,  0,  kappaw, 0];

LL = [ 0, 0, 0, -sigma,
		 0, 0, 0, 0,
		 0, 0, 0, 0];

MM = [ 0, 0, 0, sigma,
		 0, 0, kappap, 0,
		 0, 0, 0, 0];

NN = [rhoa, 0, 0, 0,
		 0,   0, 0, 0,
		 0,   0, 0, 0,
		 0,   0, 0, rhog];

% Setting the options:

[l_equ,m_states] = size(AA);
[l_equ,n_endog ] = size(CC);
[l_equ,k_exog  ] = size(DD);


DISPLAY_IMMEDIATELY=0; % Et dona els missatges d'error (en plan si no hi ha solucio al sistema, etc....
MANUAL_ROOTS=0;		  % Aixo en principi no es fa servir pero es el default
IGNORE_VV_SING = 1;    % Idem  
PROBLEM=0;				  % Variable binaria, si troba solucio PROBLEM=0; si no la troba, PROBLEM=1


TOL = .000001;         % Roots smaller than TOL are regarded as zero.
                       % Complex numbers with distance less than TOL are regarded as equal.
  
% Starting the calculations:

solve2;