function pit = varestimate(yraw,xraw,actraw,oosraw)
% Estimates a Bayesian VAR, by taking the var-covar as known a priori, and
% imposing a conditional prior on the coefficients
% Formulas are as in Giannone, Lenza, Primiceri (2011), page 2
y = []; x = []; oos = []; act = [];

% construct data series based on the variables that do have observations
% and are not an NaN
for i = 1:size(yraw,2)
    if ((sum(isfinite(yraw(:,i))==1) >= size(yraw,1)) && (sum(isfinite(actraw(:,i))==1) >= size(actraw,1)))
        y = [y yraw(:,i)];
        act = [act actraw(:,i)];
   end
end

for i = 1:size(xraw,2)
    if ((sum(isfinite(xraw(:,i))==1) >= size(xraw,1)) && (sum(isfinite(oosraw(:,i))==1) >= size(oosraw,1)))
        x = [x xraw(:,i)];
        oos = [oos oosraw(:,i)];
    end
end


[t n] = size(y);
[t k] = size(x);
lag = (k-1)/n;

% parameterizes the prior
lambda = 0.2;
sigma = diag(var(y));
b = zeros(k,n);

omega = zeros(k-1);
for l = 1:lag
    omegatemp = zeros(n);
    for i = 1:n
        for j = 1:n
            if j == i
                omegatemp(i,j) = lambda^2/l^2;
            else
                omegatemp(i,j) = lambda^2/(l^2*sigma(j,j));
            end
        end
    end
    omega(l*n-n+1:l*n,l*n-n+1:l*n) = omegatemp;
end

omega = [lambda^2     zeros(1,k-1);
         zeros(k-1,1) omega];

% calculate the posteriors
betahat = pinv(x'*x + pinv(omega))*(x'*y + pinv(omega)*b);
vhat1 = kron(sigma(1,1),pinv(x'*x + pinv(omega)));
vhat2 = kron(sigma(2,2),pinv(x'*x + pinv(omega)));

pit.pitgdp = normcdf(act(1,1),oos*betahat(:,1), oos*vhat1*oos' + sigma(1,1));
pit.pitpgdp = normcdf(act(1,2),oos*betahat(:,2), oos*vhat2*oos' + sigma(2,2));

pit.scoregdp = normpdf(act(1,1),oos*betahat(:,1), oos*vhat1*oos' + sigma(1,1));
pit.scorepgdp = normpdf(act(1,2),oos*betahat(:,2), oos*vhat2*oos' + sigma(2,2));