function [vardraw, r_pt2, q_pt2] = draw_var(x,r_pt,q_pt,r_m1,r_m2,q1,q2,r_sig,var_min,tau1,tau2)
% Construct draw of variance

n = size(x,1);
vague = 1000;
r_p=.086;
q_p=1.0;
bsum = tril(ones(n+1));
r_sig2 = r_sig^2;

% construct estimate of var_eps
lnres2 = log(x.^2);

% initial draws of Indicator Variables
tmp = rand(n,1);
ir = tmp < r_pt;
tmp = rand(n,1);
iq = tmp < q_pt;

% compute system parameters given indicators
mut = (ir*r_m1) + ((1-ir)*r_m2);
qt = (iq*q1) + ((1-q1)*q2);

% Compute Covariance Matrix
vd = diag([vague; qt]);
valpha = bsum*vd*bsum';
vy = valpha(2:n+1,2:n+1);
cy = valpha(1:n+1,2:n+1);
vy = vy + r_sig2*eye(n);
vyi = eye(n)/vy;
kgain = cy*vyi;

% Compute draws of state and shocks
ye = lnres2 - mut;
ahat0 = kgain*ye;
% ahat1 = ahat0(2:n+1);
vhat0 = valpha - kgain*cy';   % (n+1)x(n+1)
cvhat0 = chol(vhat0,'lower');
adraw0 = ahat0 + cvhat0 * randn(n+1,1);
adraw1 = adraw0(2:n+1);
edraw = lnres2 - adraw1;
udraw = adraw0(2:n+1) - adraw0(1:n);

% compute mixture probability
% e shocks ( note sigma is the same)
f1 = exp( -0.5*(((edraw - r_m1)/r_sig).^2) );
f2 = exp( -0.5*(((edraw - r_m2)/r_sig).^2) );
fe = r_p*f1 + (1-r_p)*f2;
r_pt2 = (r_p*f1)./fe;

% u shocks --- means are both zero
f1 = 1/tau1 * exp(-0.5*((udraw/tau1).^2) );
f2 = 1/tau2 * exp(-0.5*((udraw/tau2).^2) );
fu = q_p*f1 + (1-q_p)*f2;
q_pt2 = (q_p*f1)./fu;

% compute variance of AR shocks and variance of annual difs
vardraw = exp(adraw1);  % nx1

% Impose minimum value
tmp = var_min*ones(size(vardraw,1),1);
vardraw = max([vardraw tmp], [], 2);

end


