function [impulses_news,DECOMP_FEVD,NEWDECOMP_FEVD,news_series,sur_series, impulses_sur] = bvarnews(Y,p,ndraws,nhor,ident)

[T, M] = size(Y);
K = M*p+1;
nhor1=60;
% 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 = [0,0,0,0,0,0,0,0,0];
% lev = [1,1,1,1,1,0,0,0,0];
% lev = [1,1,1,1,1,1];
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.5;[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);  

DECOMP_FEVD =  zeros(ndraws,M,nhor);
NEWDECOMP_FEVD = zeros(ndraws,M,nhor);
news_series = zeros(T,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));
        irfshock = [];
    elseif ident >= 2
        if ident == 2
            V9 = IRF((M+1:M:(M-1)*M+1),1:nhor)';
            WWW = zeros(M-1);
            for j = nhor   % This corresponds to MFEVD method 
                 WWW = WWW + V9(j,:)'*V9(j,:);
            end
        elseif ident == 3
            V9 = IRF((M+1:M:(M-1)*M+1),1:nhor)';
            WWW = zeros(M-1);
            for j = 1:nhor   % This corresponds to BS method      
                WWW = WWW + (nhor+1-j)*V9(j,:)'*V9(j,:);
            end  
        elseif ident == 4
            V9 = IRF((1:M:(M-1)*M+M),1:nhor)';
            WWW = zeros(M);
            for j = nhor   % This corresponds to KS method 
                 WWW = WWW + V9(j,:)'*V9(j,:);
            end
        end
        
        if ident == 4 % for KS identification
            [Q,QQ] = eig(WWW);       % Decompose WWW
            RQ = diag(QQ);           % Main diagonal has eigenvalues
            [~,maxi] = max(RQ);   % Maximum is first eigenvalue   
            gamma = [Q(:,maxi)'];    % Take the 1st eigenvector (PC estimate)
        else % for all other identifications with zero impact restriction
            [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)
        end
               
        newshock = zeros(M,nhor);
        irfshock = zeros(M,nhor);
        
        for i = 1:M
            newshock(i,:) = IRF((0:M:(M-1)*M)+i,1:nhor)'*gamma';
            irfshock(i,:) = IRF(i,1:nhor)';
        end
        
        newshock = sign(newshock(1,nhor))*newshock;
        gamma = sign(newshock(1,nhor))*gamma;
        
    end
    
    % Calculate variance decomposition
    SIRF = cumsum(newshock.^2,2);        
    NEWDECOMP = SIRF./TOT';
    
%     % Set selection vectors for different variables in the VAR
%     pspread = zeros(1,M); pshort = zeros(1,M); pinf = zeros(1,M); pshortaug = zeros(1,M*p); pinfaug = zeros(1,M*p);
%     pspread(6) = 1; pshort(7) = 1; pinf(5) = 1; pshortaug(7) = 1; pinfaug(5)=1;
%     
%     % Long rate shock
%     LR = (pspread + pshort)*newshock;
%     RSR = (pshort - pinf)*newshock;
%     SR = pshort*newshock;
%     for ihor = 1:nhor        
%         RSR2(1,ihor) = SR(1,ihor) - pinfaug*((B^(ihor-1))*[shock*gamma ; zeros(M*p-M,1)]);
%         LEH(1,ihor) = (1/bondper)*(pshortaug*inv(eye(K-1)-B)*(eye(K-1)-B^bondper))*((B^(ihor-1))*[shock*gamma ; zeros(M*p-M,1)]);
%         RLEH(1,ihor) = LEH(1,ihor) - pinfaug*((B^(ihor-1))*[shock*gamma ; zeros(M*p-M,1)]);
%         SpreadEH(1,ihor) = LEH(1,ihor) - pshortaug*((B^(ihor-1))*[shock*gamma ; zeros(M*p-M,1)]);
%     end
%     TP = LR - LEH;
%     
%     lehfevd = cumsum(LEH.^2);
%     tpfevd = cumsum(TP.^2);
%     lrfevd = cumsum(LR.^2);
%     
%     DEC1 = tpfevd./lrfevd;
%     DEC2 = lehfevd./lrfevd;
    impulses_news(:,:,irep) = [newshock; ]; 
    impulses_sur(:,:,irep) = [irfshock; ]; 
    DECOMP_FEVD(irep,:,:) = DECOMP(1:M,1:nhor);
    if ident == 1
    news_series(:,irep) = res*gamma;
    elseif ident == 2
    news_series(:,irep) = res*gamma';
    end
    sur_series(:,irep) = res(:,1);
    NEWDECOMP_FEVD(irep,:,:) = NEWDECOMP(:,1:nhor);
end


