function [news_shock,irf_shock,news_decomp,un_decomp,x] = newshock_VAR(X,A,B,U,V,IRF,p,q,T,h,hh)

[decomp,tot] = vardecomp_eric(IRF,h,p,q);


V9 = IRF((q+1:q:(q-1)*q+1),1:hh)';

WWW = zeros(q-1);
for j = 8:hh
    WWW = WWW + (hh+1-j)*V9(j,:)'*V9(j,:);
end

[Q,QQ] = eig(WWW);
RQ = diag(QQ);
[~,maxi] = max(RQ);
x = [0 Q(:,maxi)'];

newshock = zeros(q,hh);
for i = 1:q
    newshock(i,:) = IRF((0:q:(q-1)*q)+i,1:hh)'*x';           
end

if newshock(1,hh)<0
    newshock = -1*newshock;
    x = -x;
end

% calculate variance decomposition
SIRF = newshock.^2; % This squares every element of the IRF
for i = 1:q
    SIRF(i,:) = cumsum(SIRF(i,:));
end
newsdecomp = zeros(q,hh);
for i = 1:q
    for j = 1:hh
        newsdecomp(i,j) = SIRF(i,j)/tot(j,i);
    end
end
%================| BOOTSTRAP |==================
% Compute 1 std deviation bands for IRFs
% option I
R = 5000; RR= 5000; pp=16;

[BIAS] = boot_firststage(A,U,q,p,T,V,hh,IRF,B,X,R);
[SED SEDalt low high] = boot_second(A,U,q,p,T,V,h,hh,IRF,B,BIAS,X,RR,pp);

% newshocklow = newshock - SED;
% newshockhigh = newshock + SED;

% option II
IRFlow = IRF(:,1:hh) - SEDalt;
IRFhigh = IRF(:,1:hh) + SEDalt;
[decomplow, totlow] = vardecomp_eric(IRFlow,h,p,q);
[decomphigh, tothigh] = vardecomp_eric(IRFhigh,h,p,q);

V9low = low((q+1:q:(q-1)*q+1),1:hh)' ;
V9high = high((q+1:q:(q-1)*q+1),1:hh)' ;
WWWlow = zeros(q-1);
WWWhigh = zeros(q-1);
for j = 8:hh
    WWWlow = WWWlow + (hh+1-j)*V9low(j,:)'*V9low(j,:);
    WWWhigh = WWWhigh + (hh+1-j)*V9high(j,:)'*V9high(j,:);
end

[Qlow,QQlow] = eig(WWWlow);
[Qhigh,QQhigh] = eig(WWWhigh);


RQlow = diag(QQlow);
RQhigh = diag(QQhigh);
[~,maxilow] = max(RQlow);
[~,maxihigh] = max(RQhigh);
xlow = [0 Qlow(:,maxilow)'];
xhigh = [0 Qhigh(:,maxihigh)'];

newshocklow = zeros(q,hh);
newshockhigh = zeros(q,hh);
for i = 1:q
    newshocklow(i,:) = IRFlow((0:q:(q-1)*q)+i,1:hh)'*x'; 
    newshockhigh(i,:) = IRFhigh((0:q:(q-1)*q)+i,1:hh)'*x';
end

% calculate variance decomposition
SIRFlow = newshocklow.^2; % This squares every element of the IRF
SIRFhigh = newshockhigh.^2; % This squares every element of the IRF
for i = 1:q
    SIRFlow(i,:) = cumsum(SIRFlow(i,:)); 
    SIRFhigh(i,:) = cumsum(SIRFhigh(i,:));
end
newsdecomplow = zeros(q,hh);
newsdecomphigh = zeros(q,hh);
for i = 1:q
    for j = 1:hh
        newsdecomplow(i,j) = SIRFlow(i,j)/totlow(j,i);
        newsdecomphigh(i,j) = SIRFhigh(i,j)/tothigh(j,i);
    end
end

news_shock(:,:,1) = newshocklow;
news_shock(:,:,2) = newshock;
news_shock(:,:,3) = newshockhigh;

irf_shock(:,:,1) = IRFlow;
irf_shock(:,:,2) = IRF(:,1:hh);
irf_shock(:,:,3) = IRFhigh;

news_decomp(:,:,1) = newsdecomplow;
news_decomp(:,:,2) = newsdecomp;
news_decomp(:,:,3) = newsdecomphigh;

un_decomp(:,:,1) = decomplow;
un_decomp(:,:,2) = decomp;
un_decomp(:,:,3) = decomphigh;

% total(:,:,1) = totlow;
% total(:,:,2) = tot;
% total(:,:,3) = tothigh;
