function [impulses_news] = bvarnews(Y,p,ndraws,nhor,ident,bondper)

[T, M] = size(Y);
K = M*p+1;

% Generate lagged Y matrix.
ylag = mlag2(Y,p); 
ylag = ylag(p+1:end,:);

% Redefine VAR variables y
y = Y(p+1:T,:);
% Time series observations
T = size(y,1); %#ok<*NASGU>
x = [ylag ones(T,1)];

lev = [1,0,1,1,1,1 1];
%lev = [1,1,1,1,0,0,0];
prior = 1;
lambda1 = [1e-4:1e-4:9e-4, 0.001:0.001:0.009, 0.01:0.01:0.1, 0.1:0.1:1];
lambda2 = [0.01,0.05,0.1,0.5,1,5];

lambda(:,1) = repmat(lambda1,1,length(lambda2))';
g=kron(lambda2,ones(length(lambda1),1));
lambda(:,2) = g(:);
for nmod = 1:length(lambda)
    [hm,bm] = minneprc(y,x,p,1,1,lev,prior,lambda(nmod,:)); %Minnesota prior
    xxx = inv(kron(eye(M),x'*x) + diagrv(eye(size(hm,1)),hm));  
    bb = xxx * (vec(x'*y) + diagrv(eye(size(hm,1)),hm)*bm);    %stacked posterior means of regression coeffs: vec(bpost) = inv[inv(H) + kron(I,X'X)] * [vec[X'Y] + inv(H)*bprior]
    b=reshape(bb,K,M);    %posterior coefficient matrix 
    res = y - x*b;                  %T x nvar matrix of residuals
    vmat = (1/T)*(res'*res);

    BIC(nmod) = log(det(vmat)) + (log(T)/T)*numel(b);
end
% Choose best model
[~,opt_mod] = min(BIC);

% Estimate news shock based on the optimal model
[hm,bm] = minneprc(y,x,p,1,1,lev,prior,lambda(opt_mod,:)); %Minnesota prior
xxx = inv(kron(eye(M),x'*x) + diagrv(eye(size(hm,1)),hm));  
bb = xxx * (vec(x'*y) + diagrv(eye(size(hm,1)),hm)*bm);    %stacked posterior means of regression coeffs: vec(bpost) = inv[inv(H) + kron(I,X'X)] * [vec[X'Y] + inv(H)*bprior]
b=reshape(bb,K,M);    %posterior coefficient matrix 
res = y - x*b;                  %T x nvar matrix of residuals
vmat = (1/T)*(res'*res);

sxx=chol(xxx)';
sinv=chol(inv(vmat));   
datadraws=[];
vdraws=[];
randn('state',0);
randmatrix = randn(M,T,ndraws);  
randvec = randn(K*M,ndraws);  

for irep = 1:ndraws
    if mod(irep,100)==0
        disp(['This is iteration ' num2str(irep)])
    end
       
    [bbdraw,vmatdraw] = niw(b,sxx,sinv,T,randmatrix(:,:,irep),randvec(:,irep)); %drawing from posterior coefficient matrix           
    bdraw=reshape(bbdraw,K,M);           
    res = y - x*bdraw;

    B = [bdraw(1:end-1,:)' ; eye(M*(p-1)) zeros(M*(p-1),M)];
       
    % st dev matrix for structural VAR       
    shock = chol(vmatdraw)';
       
    % 2) Estimate IRF based on Cholesky
    IRF = IRFVAR(B,shock,p,nhor);
            
    % 3) Perform variance decomposition based on IRF   
    [DECOMP,TOT] = vardecomp_eric(IRF,nhor,p,M);
    
    % 4) Extract News Shock
    if ident == 1 %Short+Long restrictions
        CC = zeros(M,M,nhor+1);
        for ii = 1:nhor+1; CC(:,:,ii) = reshape(IRF(:,ii),M,M); end;                
        SurpriseImp = CC(:,1,:);
        options = optimset('MaxFunEval',2000000,'MaxIter',1000000,'TolX',1e-15);
        [thetap, ~, flag] = fminsearch(@LongRunEffect,sign(CC(1,2:end,end))'.*ones(M-1,1),options,squeeze(CC(1,:,end)));
        gamma=[0;thetap/norm(thetap)];
        [NewsImp] = ComputeIrfOneShock(CC,gamma); %news shock
        newshock = squeeze(NewsImp(:,:,1:nhor));
        
    elseif ident == 2
        V9 = IRF((M+1:M:(M-1)*M+1),1:nhor)';
     
        WWW = zeros(M-1);
        for j = 5:nhor         
            WWW = WWW + (nhor+1-j)*V9(j,:)'*V9(j,:);
        end
        
        [Q,QQ] = eig(WWW);       % Decompose WWW
        RQ = diag(QQ);           % Main diagonal has eigenvalues
        [maxv,maxi] = max(RQ);   % Maximum is first eigenvalue   
        gamma = [0 Q(:,maxi)'];    % Take the 1st eigenvector (PC estimate)
       
        newshock = zeros(M,nhor);
        for i = 1:M
            newshock(i,:) = IRF((0:M:(M-1)*M)+i,1:nhor)'*gamma';           
        end
        
        newshock = sign(newshock(1,nhor))*newshock;
        gamma = sign(newshock(1,nhor))*gamma;                
    end
    
    
    
    impulses_news(:,:,irep) = [newshock; ];    
end


