function xtilda = DKdrawsGDPUnempStateAR1(data,param,T,n)
%% This function obtains draws of the states according to Durbin and 
% Koopman (2002)
%
% y(:,t) = dd(:,t) + ZZ(:,:,t) a(:,t) + WW(:,:,t)*eps(:,t),  eps(:,t) ~ N(0,I)
% a(:,t+1) = cc(:,t) + TT a(:,t) + RR*eta(:,t), eta(:,t) ~ N(0,I)
% a(1) ~ N(a1,C1*C1')
%
% INPUTS:
% y - data, (pxT) (where p is the number of observables and T is the
%     number of time periods)
% cc (pxT), ZZ (pxkxT), WW (pxpxT) - parameters of the observation equation
% dd (kxT), TT (kxk), RR (kxk) - parameters of the state equation
% a1 (kx1), C1 (kxk) - parameters of the distribution of the initial state
% p (scalar) - number of observables
% k (scalar) - number of states
% T (scalar) - sample size
% ind_output - 1: return the smoothed states and the log likelihood [aaa,loglik]
%              2: return a draw of the states from the simulation smoother and zeros as the log likelihood [aaa,zeros(1,T)]
% OUTPUTS:
% if ind_output = 1
%   aaa = the mean of the states conditional on y, (kxT)
%   loglik = log likelihood of each observation of y, (1xT)
% if ind_output = 2
%   aaa = a draw of the states conditional on y, (kxT)
%   loglik = []
%
% Manuel Gonzalez-Astudillo
% 11/16/2017

%% Parameterizing
phi1 = param(1);
phi2 = param(2);
mui = param(3:n+3-1);
alphai = param(n+3:2*n+3-1);
deltayi = param(2*n+3:3*n+3-1);
deltaui = param(3*n+3:4*n+3-1);
theta1i = param(4*n+3:5*n+3-1);
theta2i = param(5*n+3:6*n+3-1);
rhoi = param(6*n+3:7*n+3-1);
sigmazetai2 = param(7*n+3:8*n+3-1);
sigmaetayi2 = param(8*n+3:9*n+3-1);
sigmaetaui2 = param(9*n+3:10*n+3-1);
sigmae2 = 1;
sigmaetay2 = 1;
sigmaetau2 = 1;
mu = 0;

%% Setting up data
y = data;

%% Setting up matrices
% The order of the latent variables is:
% c_t, c_{t-1}, c_{t-2}, eta_{yt}, eta_{ut}, nu_{it}, nu_{it-1}, nu_{it-2}
k1 = 5; % There are 5 common components
ks1 = 0;
ks2 = 3;
ks = ks1 + ks2; % There are 3 latent variables nu_it (t,t-1,t-2)
k2 = ks*n;
k = k1+k2; % There are 5 common factors and 3 latent variables for each of the 51 states

ZZdUsmall = zeros(n,k2);
for i = 1:n
    ZZdUsmall(i,(i-1)*ks+1:i*ks) = [theta1i(i) (theta2i(i)-theta1i(i)) -theta2i(i)];
end
ZZ = [alphai -alphai zeros(n,1) deltayi zeros(n,1) kron(eye(n),[1 -1 0]);
      theta1i.*alphai alphai.*(theta2i-theta1i) -alphai.*theta2i  zeros(n,1) deltaui ZZdUsmall];
dd = [mui+mu*deltayi;zeros(n,1)];
HH = diag([sigmaetayi2; sigmaetaui2]);
WW = sqrt(HH);

cc = zeros(k,1);
TT_11 = [phi1 phi2 0 0 0;
         1 0 0 0 0;
         0 1 0 0 0;
         0 0 0 0 0;
         0 0 0 0 0];
TT_12 = zeros(k1,k2);
TT_21 = zeros(k2,k1);
TT_22small = [0 0 0;
              1 0 0;
              0 1 0];
TT_22 = kron(eye(n),TT_22small);

TT = [TT_11 TT_12; TT_21 TT_22];

Rhodiag = zeros(k,1);
Rhodiag(k1+ks1+1:ks:end) = rhoi;
Rho = diag(Rhodiag);
TT = TT + Rho;

RQRdiag = zeros(k,1);
RQRdiag(1) = sigmae2;
RQRdiag(4) = sigmaetay2;
RQRdiag(5) = sigmaetau2;
RQRdiag(k1+ks1+1:ks:end) = sigmazetai2;
RQR = diag(RQRdiag);
RR = sqrt(RQR);

a1 = zeros(k,1);
P1diag = zeros(k,1);
P1diag(1:3) = sigmae2*(1-phi2)/((1+phi2)*((1-phi2)^2-phi1^2));
P1diag(4) = sigmaetay2;
P1diag(5) = sigmaetau2;
P1diag(k1+ks1+1:ks:end) = sigmazetai2./(1-rhoi.^2);
P1diag(k1+ks1+2:ks:end) = sigmazetai2./(1-rhoi.^2);
P1diag(k1+ks1+3:ks:end) = sigmazetai2./(1-rhoi.^2);
C1 = sqrt(diag(P1diag));

tVec = ones(T,1);
dd = dd(:,tVec);
ZZ = ZZ(:,:,tVec);
WW = WW(:,:,tVec);
cc = cc(:,tVec);

xtilda = dksimulator(y,dd,ZZ,WW,cc,TT,RR,a1,C1,2);

end

