global endo_nbr ykmin_ xkmin_ it_ options_ iy_
% nx = size(dr.ghu,2);
% k1=dr_.order_var;
% A=dr_.ghx(k1,k1(klag([1:nx],1)));
% B=dr_.ghu(k1,:);
% ys=dr_.ys;

gammax=oo_.SmoothedShocks.gammax;
zetacx=oo_.SmoothedShocks.zetacx;
zetaix=oo_.SmoothedShocks.zetaix;
PATRIOT=oo_.SmoothedShocks.PATRIOT;
mps=oo_.SmoothedShocks.mps;
lambdafx=oo_.SmoothedShocks.lambdafx;

datagammaU=oo_.SmoothedVariables.gammaU+gammaUU;
dataepsilU=oo_.SmoothedVariables.epsilU+epsilUU;
datazetacU=oo_.SmoothedVariables.zetacU+zetacUU;
datazetaiU=oo_.SmoothedVariables.zetaiU+zetaiUU;
datalambdafU=oo_.SmoothedVariables.lambdafU+lambdafUU;

v= [gammax(1:end-1) lambdafx(1:end-1) mps(1:end-1) PATRIOT(1:end-1) zetacx(1:end-1) zetaix(1:end-1)];
CC=corrcoef(v);

zrs=zeros(1,length(gammax)-1);

is=11;%if between 1 and 10, simulate response to just isth shock...if > 10 do all

iss=[1 2 3 6 9 10];

for ij = 1:length(iss)
    
    u=zeros(10,length(gammax)-1);
    
    is=iss(ij);

    figure('name',['shock ' num2str(is)]);

    if is == 1|is > 10;u(1,:)=gammax(1:end-1)';end
    if is == 2|is > 10;u(2,:)=lambdafx(1:end-1)';end
    if is == 3|is > 10;u(3,:)=mps(1:end-1)';end
    if is == 6|is > 10;u(6,:)=PATRIOT(1:end-1)';end
    if is == 9|is > 10;u(9,:)=zetacx(1:end-1)';end
    if is == 10|is > 10;u(10,:)=zetaix(1:end-1)';end

    structnames=fieldnames(oo_.SmoothedVariables);
    yt=zeros(size(structnames,1),1);
    for ii = 1:size(structnames,1)
        yt(ii)=oo_.SmoothedVariables.(char(structnames(ii,:)))(1,1);
    end

    %the variables in yt are in deviation from steady state and ordered....
    yt=yt(dr_.inv_order_var)+ys_;

    iorder=options_.order;
    ex=u';

    yy = simult_(yt,dr_,ex,iorder);

    tt=[ykmin_:211+ykmin_-1];
    
    readdatalhu

    gammaU=yy(8,:);
    epsilU=yy(5,:);
    zetacU=yy(40,:);
    zetaiU=yy(41,:);
    lambdafU=yy(11,:);

    subplot(431)
    plot(time,cumsum(yy(4,tt)),time,cumsum(dly),time,cumsum(ones(length(dly),1)*log(muzstarUU)))
    legend('simulation','data')
    title('dly')
    axis tight
    subplot(432)
    plot(time,yy(13,tt),time,lhU,time,ones(length(time),1)*log(hUU))
    title('lhours')
    axis tight
    subplot(433)
    plot(time,yy(14,tt),time,liy,time,ones(length(time),1)*log(iUU/(iUU+cUU)))
    title('liy')
    axis tight
    subplot(434)
    plot(time,cumsum(yy(3,tt)),time,cumsum(dlDJ),time,cumsum(ones(length(dly),1)*log(muzstarUU)))
    title('dlDJ')
    axis tight
    subplot(435)
    plot(time,yy(29,tt),time,RXU,time,ones(length(dly),1)*RXUU)
    title('interest rate')
    axis tight
    subplot(436)
    plot(time,yy(22,tt),time,piU,time,ones(length(dly),1)*piUU)
    title('inflation')
    axis tight
    subplot(437)
    plot(time,gammaU(tt),time,datagammaU,time,ones(length(tt),1)*gammaUU)
    title('gamma')
    axis tight
    subplot(438)
    plot(time,epsilU(tt),time,dataepsilU)
    title('epsil')
    axis tight
    subplot(439)
    plot(time,zetacU(tt),time,datazetacU)
    title('zetac')
    axis tight
    subplot(4,3,10)
    plot(time,zetaiU(tt),time,datazetaiU)
    title('zetai')
    axis tight
    subplot(4,3,11)
    plot(time,lambdafU(tt),time,datalambdafU)
    title('lambdaf')
    axis tight
end

figure

analshocks