function [eps, eta, tau] = draw_eps_eta(y, var_eps_n, var_eta_n)
% Construct draws of eps amd eta

n = size(y,1);
ltone = tril(ones(n));
cov_eps = diag(var_eps_n);
cov_tau = ltone*cov_eps*ltone';
cov_y = cov_tau + diag(var_eta_n);
kappa = cov_tau/cov_y;
mutau_y = kappa*y;
covtau_y = cov_tau - kappa*cov_tau';
chol_covtau_y = chol(covtau_y,'lower');
tau = mutau_y + chol_covtau_y * randn(n,1);
eta = y - tau;
eps = [tau(1); diff(tau)];
end

