
function z=usmodel_pshk_fff(y)
z=zeros(41,1);
global ex_ ex_det_ it_ recur_

global calfa cbeta cbetabar ccs ccy cfc cg cgamma cgy chabb 
global cik cikbar cindp cindw cinvs ciy cky clandap clandaw clk 
global cmap cmaw constebeta constelab constepinf conster cpie cprobp cprobw ...
 cr 
global crdpi crdy crhoa crhoas crhob crhog crhols crhoms crhopinf crhoqs 
global crhow crk crkky crpi crr cry csadjcost csigl csigma ctou 
global ctrend curvp curvw cw cwhlc cwly czcap 
usmodel_stst;
z(1) = 0 -(calfa*y(31)+(1-calfa)*(y(37))-0*(1-calfa)*y(1)-1*y(1)+1/((1/(1+ ...
cbetabar*cgamma*cindp))*((1-cprobp)*(1-cbetabar*cgamma*cprobp)/cprobp)/((cfc ...
-1)*curvp+1))*y(34));
z(2) = y(41) -((1/(czcap/(1-czcap)))*y(31));
z(3) = y(31) -((y(37))+y(19)-y(15));
z(4) = y(15) -(y(17)+y(41));
z(5) = y(13) -((1/(1+cbetabar*cgamma))*(y(13)+cbetabar*cgamma*y(13)+(1/( ...
cgamma^2*csadjcost))*y(27))+y(28));
z(6) = y(27) -(-y(33)-0*y(2)+(1/((1-chabb/cgamma)/(csigma*(1+chabb/cgamma)) ...
))*y(2)+(crk/(crk+(1-ctou)))*y(31)+((1-ctou)/(crk+(1-ctou)))*y(27));
z(7) = y(4) -((chabb/cgamma)/(1+chabb/cgamma)*y(4)+(1/(1+chabb/cgamma))* ...
y(4)+((csigma-1)*cwhlc/(csigma*(1+chabb/cgamma)))*(y(19)-y(19))-(1-chabb/ ...
cgamma)/(csigma*(1+chabb/cgamma))*(y(33)+0*y(2))+y(2));
z(8) = y(39) -(ccy*y(4)+ciy*y(13)+y(11)+crkky*y(41));
z(9) = y(39) -(cfc*(calfa*y(15)+(1-calfa)*y(19)+y(1)));
z(10) = y(37) -(csigl*y(19)+(1/(1-chabb/cgamma))*y(4)-(chabb/cgamma)/(1- ...
chabb/cgamma)*y(4)+1/((1-cprobw)*(1-cbetabar*cgamma*cprobw)/((1+cbetabar* ...
cgamma)*cprobw)*(1/((clandaw-1)*curvw+1)))*y(21));
z(11) = y(17) -((1-cikbar)*y(17)+(cikbar)*y(13)+(cikbar)*(cgamma^2* ...
csadjcost)*y(28));
z(12) = y(22) -(calfa*y(30)+(1-calfa)*(y(36))-1*y(1)-0*(1-calfa)*y(1));
z(13) = y(40) -((1/(czcap/(1-czcap)))*y(30));
z(14) = y(30) -(y(36)+y(18)-y(14));
z(15) = y(14) -(y(16)+y(40));
z(16) = y(12) -((1/(1+cbetabar*cgamma))*(y(12)+cbetabar*cgamma*y(12)+(1/( ...
cgamma^2*csadjcost))*y(26))+y(28));
z(17) = y(26) -(-y(29)+y(24)-0*y(2)+(1/((1-chabb/cgamma)/(csigma*(1+chabb/ ...
cgamma))))*y(2)+(crk/(crk+(1-ctou)))*y(30)+((1-ctou)/(crk+(1-ctou)))*y(26));
z(18) = y(3) -((chabb/cgamma)/(1+chabb/cgamma)*y(3)+(1/(1+chabb/cgamma))* ...
y(3)+((csigma-1)*cwhlc/(csigma*(1+chabb/cgamma)))*(y(18)-y(18))-(1-chabb/ ...
cgamma)/(csigma*(1+chabb/cgamma))*(y(29)-y(24)+0*y(2))+y(2));
z(19) = y(38) -(ccy*y(3)+ciy*y(12)+y(11)+1*crkky*y(40));
z(20) = y(38) -(cfc*(calfa*y(14)+(1-calfa)*y(18)+y(1)));
z(21) = y(24) -((1/(1+cbetabar*cgamma*cindp))*(cbetabar*cgamma*y(24)+cindp* ...
y(24)+((1-cprobp)*(1-cbetabar*cgamma*cprobp)/cprobp)/((cfc-1)*curvp+1)*( ...
y(22)))+y(34));
z(22) = y(36) -((1/(1+cbetabar*cgamma))*y(36)+(cbetabar*cgamma/(1+cbetabar* ...
cgamma))*y(36)+(cindw/(1+cbetabar*cgamma))*y(24)-(1+cbetabar*cgamma*cindw)/( ...
1+cbetabar*cgamma)*y(24)+(cbetabar*cgamma)/(1+cbetabar*cgamma)*y(24)+(1- ...
cprobw)*(1-cbetabar*cgamma*cprobw)/((1+cbetabar*cgamma)*cprobw)*(1/((clandaw ...
-1)*curvw+1))*(csigl*y(18)+(1/(1-chabb/cgamma))*y(3)-((chabb/cgamma)/(1- ...
chabb/cgamma))*y(3)-y(36))+1*y(35));
z(23) = y(29) -(crpi*(1-crr)*y(24)+cry*(1-crr)*(y(38)-y(39))+crdy*(y(38)- ...
y(39)-y(38)+y(39))+crr*y(29)+y(23));
z(24) = y(1) -(crhoa*y(1)+ex_(it_-1,1));
z(25) = y(2) -(crhob*y(2)+ex_(it_-1,2));
z(26) = y(11) -(crhog*(y(11))+ex_(it_-1,3)+cgy*ex_(it_-1,1));
z(27) = y(28) -(crhoqs*y(28)+ex_(it_-1,6));
z(28) = y(23) -(crhoms*y(23)+ex_(it_-1,4));
z(29) = y(34) -(crhopinf*y(34)+y(9)-cmap*y(9));
z(30) = y(9) -(ex_(it_-1,5));
z(31) = y(21) -(crhow*y(21)+ex_(it_-1,7));
z(32) = y(35) -(y(21)+ex_(it_-1,8));
z(33) = y(10) -(0);
z(34) = y(16) -((1-cikbar)*y(16)+cikbar*y(12)+cikbar*cgamma^2*csadjcost* ...
y(28));
z(35) = y(8) -(y(38)-y(38)+ctrend);
z(36) = y(5) -(y(3)-y(3)+ctrend);
z(37) = y(6) -(y(12)-y(12)+ctrend);
z(38) = y(7) -(y(36)-y(36)+ctrend);
z(39) = y(25) -(1*(y(24))+constepinf);
z(40) = y(32) -(1*(y(29))+conster);
z(41) = y(20) -(y(18)+constelab);
if ~isreal(z)
  z = real(z)+imag(z).^2;
end
