function mdata = drawq_ur(param,Sbar,P,Q,ldata);
%DRAWQ_UR   generate time series for model economy
%           mdata=drawq_ur(param,Sbar,P,Q) generates productivity
%           and hours per capita for the model economy with parameters
%           in vector param and stochastic process summarized by
%           Sbar,P,Q.

%           Ellen McGrattan, 3-20-04
%           Revised, 9-5-06

%---------------------------------------------------------------------
% 1. Default parameters for 
%    gn, gz, beta, delta, psi, sigma, theta, ghat, a_0, a, b, D, R
%
gn      = param(1);
gz      = param(2);
beta    = param(3);
delta   = param(4);
psi     = param(5);
sigma   = param(6);
theta   = param(7);
P0      = (eye(4)-P)*Sbar;

%---------------------------------------------------------------------
% 2.  Compute equilibrium
%---------------------------------------------------------------------

%
% 2a. Steady state:
%
tem        = (eye(4)-P)\P0;
zs         = exp(tem(1));
tauls      = tem(2);
tauxs      = tem(3);
gs         = exp(tem(4));
kls        = ((1+tauxs)*(1-beta*zs^(-sigma)*(1-delta))/ ...
              (beta*zs^(-sigma)*theta))^(1/(theta-1))*zs;
A          = kls^(theta-1)*zs^(-theta)-(1+gn)+(1-delta)/zs;
B          = (1-tauls)*(1-theta)*kls^theta*zs^(-theta)/psi;
ks         = (B+gs)/(A+B/kls);
cs         = A*ks-gs;
ls         = ks/kls;
ys         = (ks/zs)^theta*ls^(1-theta);
xs         = ys-cs-gs;
X0s        = [log(ks);log(zs);tauls;tauxs;log(gs);1];
X0         = [X0s;X0s];
Y0         = [log(zs);log(zs);log(ls);log(zs)];

%
% 2b. Call subroutine with residuals:
%

Z          = [log(ks);log(ks);log(ks);log(zs);log(zs);tauls;tauls;
              tauxs;tauxs;log(gs);log(gs)];
del        = max(abs(Z)*1e-5,1e-8);
for i=1:11;
  Zp       = Z;
  Zm       = Z;
  Zp(i)    = Z(i)+del(i);
  Zm(i)    = Z(i)-del(i);
  dR(i,1)  = (res_wedge_ur(Zp,param)-res_wedge_ur(Zm,param))/(2*del(i));
end;

%
% 2c. Solution:  log k[t+1] = gamma0 + gammak* log k[t] + gamma* S[t]
%

a0         = dR(1);
a1         = dR(2);
a2         = dR(3);
b0         = dR(4:2:11)';
b1         = dR(5:2:11)';
tem        = roots([a0,a1,a2]);
gammak     = tem(find(abs(tem)<1));
gamma      = -((a0*gammak+a1)*eye(4)+a0*P')\(b0*P+b1)';
gamma0     = (1-gammak)*log(ks)-gamma'*[log(zs);tauls;tauxs;log(gs)];
Gamma      = [gammak;gamma;gamma0];

%
% 2d. State-space system:   X[t+1] = A X[t] + B eps[t+1]
%                           Y[t]   = C X[t] + ome[t]
%                           ome[t] = D ome[t-1] + eta[t],  E eta eta' ~ N(0,R)
%


philh      =-(psi*ys*(1-theta)+(1-theta)*(1-tauls)*ys*(1-ls)/ls*theta+ ...
             (1-theta)*(1-tauls)*ys);
philk      = (psi*ys*theta+psi*(1-delta)*ks/zs- ...
             (1-theta)*(1-tauls)*ys*(1-ls)/ls *theta)/philh;
philz      = (psi*ys*(-theta)-psi*(1-delta)*ks/zs+ ...
             (1-theta)*theta*(1-tauls)*ys*(1-ls)/ls)/philh;
phill      = ((1-theta)*(1-tauls)*ys*(1-ls)/ls *(1/(1-tauls)))/philh;
philg      = (-psi*gs)/philh;
philkp     = (-psi*(1+gn)*ks)/philh;
phiyk      = theta+(1-theta)*philk;
phiyz      = (1-theta)*philz-theta;
phiyl      = (1-theta)*phill;
phiyg      = (1-theta)*philg;
phiykp     = (1-theta)*philkp;
phixk      = -ks/(zs*xs)*(1-delta);
phixkp     = ks/xs*(1+gn);
phixz      = ks/(zs*xs)*(1-delta);

As    = [ gammak,                         gamma', gamma0;
         [0;0;0;0],            P,                     P0;
               0,      0,      0,      0,      0,      1];
Bs    = [ 0,0,0,0;
            Q;
          0,0,0,0];
A     = [As,zeros(6);
         eye(6),zeros(6)];
B     = [Bs;
         zeros(6,4)];

phiy0 = log(ys)-phiyk*log(ks)-phiyz*log(zs)-phiyl*tauls-phiyg*log(gs)- ...
                phiykp*log(ks);
phix0 = log(xs)-phixk*log(ks)-phixz*log(zs)-phixkp*log(ks);
phil0 = log(ls)-philk*log(ks)-philz*log(zs)-phill*tauls-philg*log(gs)- ...
                philkp*log(ks);
C     = [ [phiyk-phiykp,1+phiyz,phiyl,     0,phiyg,phiy0, ...
           -phiyk,-phiyz,-phiyl,    0,-phiyg,-phiy0]+ ....
            phiykp*[Gamma(1:6)',zeros(1,6)];
          [phixk-phixkp,1+phixz,    0,     0,    0,phix0, ...
           -phixk,-phixz,     0,    0,     0,-phix0]+ ....
            phixkp*[Gamma(1:6)',zeros(1,6)];
          [philk,philz,phill,    0,philg,phil0,0,0,0,0,0,0]+ ....
            philkp*[Gamma(1:6)',zeros(1,6)];
          0,1,0,0,1,0,0,0,0,0,-1,0];
shks    = randn(ldata,4);
[Yt,Xt] = dlsim(A,B,C,zeros(4),shks,X0);
mdata   = [ys*cumprod(exp(Yt(:,1)))./exp(Yt(:,3)),exp(Yt(:,3))];

