%% Use this file for Impulse Responses

%model1 is the baseline model
%model2 is the same model with the same stochastic runs, except for an
%additional one time  shock  at a certain time point 
%run both for 500 (number of runs)
%take the difference!



function  result = impulseresponses(T,number_of_runs_irf, seed,varargin)
rng(seed,'twister');        % Set random seed
rndseed = rng;
rng(rndseed);

%% Parameters of the model
Tq=64; %number of days per quarter
h=1/Tq;%frequency 
sig=1.4; %risk aversion coefficient
c1=0.9; %stock price effect in IS equation   
rbar=0.01; %real ineterst rate in ss
betd=1/(1+h*rbar); %discount factor daily frequency
et=1; %inverse Frisch elascticity of labor supply
thet=0.55; %quarterly calvo price rigidity
thetd=1-h*(1-thet); %daily calvo price rigidity
gamd=(sig+et)*(1-thetd)*(1-betd*thetd)/thetd; 
kap=0.001; %stock price cost effect in phillips curve
alph=0.2; %extrapolation parameter in heuristics 
phi=10; %intensity of choice
for i = 1:length(varargin)
                                     if strcmpi(varargin{i},'phi')
                                     phi = varargin{i+1};
                                     end
 end
zet=0.5; %memory coefficient real
kap_s=0.04; %extrapolation for chartists and reversed speed for fundamentalists in financial sector
d=0.95; %memory financial sector
e=30; %intensity of choice financial sector
for i = 1:length(varargin)
                                     if strcmpi(varargin{i},'e')
                                     e = varargin{i+1};
                                     end
 end
a=1; %excess demand impact on stock price
k=0.2; %impact hh stock demand on excess demand
c2=0.5 ; %sensitivity to real stock price in hh stock demand
c3=1; %sensitivity to nomninal interest rate in hh stock demand
g=0.5; %bias fundamental stock price to output
for i = 1:length(varargin)
                                     if strcmpi(varargin{i},'g')
                                     g = varargin{i+1};
                                     end
 end
rho_s=0.7;  %daily stock shock persistence 
delta_x=0.5; %taylor rule output
delta_pi= 1.5; %taylor rule inflation
 for i = 1:length(varargin)
                                     if strcmpi(varargin{i},'delta_pi')
                                     delta_pi = varargin{i+1};
                                     end
 end
delta_s=0; %taylor rule stock
 for i = 1:length(varargin)
                                     if strcmpi(varargin{i},'delta_s')
                                     delta_s = varargin{i+1};
                                     end
 end
sig_s=0.01; %stock price shock standard deviation

for i = 1:length(varargin)
                                     if strcmpi(varargin{i},'sig_s')
                                     sig_s = varargin{i+1};
                                     end
end

IMPULSE_Shock=0; %extra shock  to calculate impulse responses
                                  for i = 1:length(varargin)
                                     if strcmpi(varargin{i},'IMPULSE_Shock')
                                     IMPULSE_Shock = varargin{i+1};
                                     end
                                  end
%% matrices
A=[1 0 c1+h/sig*delta_s;-gamd h kap;0 0 1+a*k*c2+a*k*c3*delta_s];
B=[1 h*(1/sig-c1) c1; 0 betd*h 0; 0 0 0];

%% start values model 1
x=zeros(T, number_of_runs_irf); %output daily
pi=zeros(T, number_of_runs_irf); %inflation daily
s=zeros(T, number_of_runs_irf); %stock price
sf=zeros(T, number_of_runs_irf); %fundamental stock price


E_pi=zeros(T, number_of_runs_irf); % expectation inflation
E_s=zeros(T, number_of_runs_irf); % expectation stock price
Et_x= zeros(T, number_of_runs_irf); % expectation output daily
Et_s_c= zeros(T, number_of_runs_irf); % expectation stock chartists
Et_s_F= zeros(T, number_of_runs_irf); % expectation stock fundamentalists
Eq_tar_x=zeros(T, number_of_runs_irf); % expectation output targetrs
Eq_tar_pi=zeros(T, number_of_runs_irf); % expectation inflation targetrs
Eq_sta_x=zeros(T, number_of_runs_irf); % expectation output static
Eq_sta_pi=zeros(T, number_of_runs_irf); % expectation inflation static
Eq_ext_x=zeros(T, number_of_runs_irf); % expectation output extrapolators
Eq_ext_pi=zeros(T, number_of_runs_irf); % expectation inflation extrapolators
Eq_x=zeros(T, number_of_runs_irf); % expectation output 
Eq_pi=zeros(T, number_of_runs_irf); % expectation inflation 


eps_s=zeros(T, number_of_runs_irf); % shock stock price

i=zeros(T, number_of_runs_irf); %interest rate daily

A_tar_pi=zeros(T, number_of_runs_irf); %attractivity targeters inflation
A_sta_pi=zeros(T, number_of_runs_irf); %attractivity static inflation
A_tar_x=zeros(T, number_of_runs_irf); %attractivity targeters output
A_sta_x=zeros(T, number_of_runs_irf); %attractivity static output
A_ext_x=zeros(T, number_of_runs_irf); %attractivity extrapolator output
A_ext_pi=zeros(T, number_of_runs_irf); %attractivity extrapolator inflation

A_f_s=zeros(T, number_of_runs_irf); %attractivity funademntalists stock price
A_c_s=zeros(T, number_of_runs_irf); %attractivity charters stock price
A_nt_s=zeros(T, number_of_runs_irf); %attractivity non-traders stock price

xq= zeros(T, number_of_runs_irf); %output quarterly value
piq= zeros(T, number_of_runs_irf); %inflation quarterly value
pq= zeros(T, number_of_runs_irf); %price level quarterly

D=zeros(T, number_of_runs_irf); % HH demand stock
D_f=zeros(T, number_of_runs_irf); % fundamentalists demand stock
D_c=zeros(T, number_of_runs_irf); % chartists demand stock


omega_x_tar=zeros (T,number_of_runs_irf); %weights output targeters
omega_pi_tar=zeros(T,number_of_runs_irf); %weights inflation targeters
omega_x_ext=zeros(T, number_of_runs_irf); %weights output extrapolators
omega_pi_ext=zeros(T, number_of_runs_irf); %weights inflation extrapolators
omega_x_sta=zeros(T, number_of_runs_irf); %weights output static
omega_pi_sta=zeros(T, number_of_runs_irf); %weights inflation static

omega_s_c=zeros(T, number_of_runs_irf); %weights stock price charters
omega_s_f=zeros(T, number_of_runs_irf); %weights stock price  fundamentalists
omega_s_nt=zeros(T, number_of_runs_irf); %weights stock price  non-traders

xqq=zeros(T,number_of_runs_irf); %output quarterly freq
piqq=zeros(T, number_of_runs_irf); %inflation quarterly freq

%% special starting values model 1
omega_x_tar(1:128,:)=1/3; 
omega_pi_tar(1:128,:)=1/3;
omega_x_sta(1:128,:)=1/3;
omega_pi_sta(1:128,:)=1/3;
omega_pi_ext(1:128,:)=1/3;
omega_x_ext(1:128,:)=1/3;
omega_s_c(1:128,:)=1/3;
omega_s_f(1:128,:)=1/3;
omega_s_nt(1:128,:)=1/3; 
q(1:64,:) = 1;
q(65:128,:) = 2;

%% start values model 2
x2=zeros(T, number_of_runs_irf); %output daily
pi2=zeros(T, number_of_runs_irf); %inflation daily
s2=zeros(T, number_of_runs_irf); %stock price
sf2=zeros(T, number_of_runs_irf); %fundamental stock price


E_pi2=zeros(T, number_of_runs_irf); % expectation inflation
E_s2=zeros(T, number_of_runs_irf); % expectation stock price
Et_x2= zeros(T, number_of_runs_irf); % expectation output daily
Et_s_c2= zeros(T, number_of_runs_irf); % expectation stock chartists
Et_s_F2= zeros(T, number_of_runs_irf); % expectation stock fundamentalists
Eq_tar_x2=zeros(T, number_of_runs_irf); % expectation output targetrs
Eq_tar_pi2=zeros(T, number_of_runs_irf); % expectation inflation targetrs
Eq_sta_x2=zeros(T, number_of_runs_irf); % expectation output static
Eq_sta_pi2=zeros(T, number_of_runs_irf); % expectation inflation static
Eq_ext_x2=zeros(T, number_of_runs_irf); % expectation output extrapolators
Eq_ext_pi2=zeros(T, number_of_runs_irf); % expectation inflation extrapolators
Eq_x2=zeros(T, number_of_runs_irf); % expectation output 
Eq_pi2=zeros(T, number_of_runs_irf); % expectation inflation 


eps_s2=zeros(T, number_of_runs_irf); % shock stock price

i2=zeros(T, number_of_runs_irf); %interest rate daily

A_tar_pi2=zeros(T, number_of_runs_irf); %attractivity targeters inflation
A_sta_pi2=zeros(T, number_of_runs_irf); %attractivity static inflation
A_tar_x2=zeros(T, number_of_runs_irf); %attractivity targeters output
A_sta_x2=zeros(T, number_of_runs_irf); %attractivity static output
A_ext_x2=zeros(T, number_of_runs_irf); %attractivity extrapolator output
A_ext_pi2=zeros(T, number_of_runs_irf); %attractivity extrapolator inflation

A_f_s2=zeros(T, number_of_runs_irf); %attractivity funademntalists stock price
A_c_s2=zeros(T, number_of_runs_irf); %attractivity charters stock price
A_nt_s2=zeros(T, number_of_runs_irf); %attractivity non-traders stock price

xq2= zeros(T, number_of_runs_irf); %output quarterly value
piq2= zeros(T, number_of_runs_irf); %inflation quarterly value
pq2= zeros(T, number_of_runs_irf); %price level quarterly

D2=zeros(T, number_of_runs_irf); % HH demand stock
D_f2=zeros(T, number_of_runs_irf); % fundamentalists demand stock
D_c2=zeros(T, number_of_runs_irf); % chartists demand stock


omega_x_tar2=zeros (T,number_of_runs_irf); %weights output targeters
omega_pi_tar2=zeros(T,number_of_runs_irf); %weights inflation targeters
omega_x_ext2=zeros(T, number_of_runs_irf); %weights output extrapolators
omega_pi_ext2=zeros(T, number_of_runs_irf); %weights inflation extrapolators
omega_x_sta2=zeros(T, number_of_runs_irf); %weights output static
omega_pi_sta2=zeros(T, number_of_runs_irf); %weights inflation static

omega_s_c2=zeros(T, number_of_runs_irf); %weights stock price charters
omega_s_f2=zeros(T, number_of_runs_irf); %weights stock price  fundamentalists
omega_s_nt2=zeros(T, number_of_runs_irf); %weights stock price  non-traders

xqq2=zeros(T,number_of_runs_irf); %output quarterly freq
piqq2=zeros(T, number_of_runs_irf); %inflation quarterly freq

%% special starting values model 2
omega_x_tar2(1:128,:)=1/3; 
omega_pi_tar2(1:128,:)=1/3;
omega_x_sta2(1:128,:)=1/3;
omega_pi_sta2(1:128,:)=1/3;
omega_pi_ext2(1:128,:)=1/3;
omega_x_ext2(1:128,:)=1/3;
omega_s_c2(1:128,:)=1/3;
omega_s_f2(1:128,:)=1/3;
omega_s_nt2(1:128,:)=1/3; 

%% shocks;
Randomz=randn(T, number_of_runs_irf);
IMPULSE_Shock_var=zeros(T, number_of_runs_irf);

%% Start dynamics!
    for MC=1:number_of_runs_irf
    
    for t=129:T
        
        %% (((shocks)))
eps_s(t,MC)  = rho_s*eps_s(t-1,MC) +sig_s*Randomz(t,MC); 


    
        
 %% fundamental value
 sf(t,MC)=g*xq(t-1,MC);
 
 
 %% Attractivity
A_tar_pi(t,MC) = zet*A_tar_pi(t-1,MC) - (Eq_tar_pi (t-2,MC)-piq(t-1,MC))^2; 
A_sta_pi(t,MC) = zet*A_sta_pi(t-1,MC) - (Eq_sta_pi (t-2,MC)-piq(t-1,MC))^2; 
A_tar_x(t,MC) = zet*A_tar_x(t-1,MC) - (Eq_tar_x (t-2,MC)-xq(t-1,MC))^2; 
A_sta_x(t,MC) = zet*A_sta_x(t-1,MC) - (Eq_sta_x (t-2,MC)-xq(t-1,MC))^2; 
A_ext_x(t,MC) = zet*A_ext_x(t-1,MC) - (Eq_ext_x (t-2,MC)-xq(t-1,MC))^2; 
A_ext_pi(t,MC) = zet*A_ext_pi(t-1,MC) - (Eq_ext_pi (t-2,MC)-piq(t-1,MC))^2; 

A_f_s(t,MC) = d*A_f_s(t-1,MC) - (exp(s (t-2,MC))-exp(s(t-1,MC)))*D_f(t-2,MC); 
A_c_s(t,MC) = d*A_c_s(t-1,MC) - (exp(s (t-2,MC))-exp(s(t-1,MC)))*D_c(t-2,MC); 
A_nt_s(t,MC)=0;

%% (((weights)))
omega_x_tar(t,MC) = exp(phi*A_tar_x(t,MC))/(exp(phi * A_tar_x(t,MC)) + exp(phi * A_ext_x(t,MC))+ exp(phi * A_sta_x(t,MC)));
omega_pi_tar(t,MC)  = exp(phi*A_tar_pi(t,MC))/(exp(phi * A_tar_pi(t,MC)) + exp(phi * A_ext_pi(t,MC))+ exp(phi * A_sta_pi(t,MC)));
omega_x_ext(t,MC)  = exp(phi*A_ext_x(t,MC))/(exp(phi * A_tar_x(t,MC)) + exp(phi * A_ext_x(t,MC))+ exp(phi * A_sta_x(t,MC)));
omega_pi_ext(t,MC)=exp(phi*A_ext_pi(t,MC))/(exp(phi * A_tar_pi(t,MC)) + exp(phi * A_ext_pi(t,MC))+ exp(phi * A_sta_pi(t,MC)));
omega_x_sta(t,MC)=exp(phi*A_sta_x(t,MC))/(exp(phi * A_tar_x(t,MC)) + exp(phi * A_ext_x(t,MC))+ exp(phi * A_sta_x(t,MC)));
omega_pi_sta(t,MC)=exp(phi*A_sta_pi(t,MC))/(exp(phi * A_tar_pi(t,MC)) + exp(phi * A_ext_pi(t,MC))+ exp(phi * A_sta_pi(t,MC)));

omega_s_nt(t,MC)=exp(e*A_f_s(t,MC))/(exp(e * A_f_s(t,MC)) + exp(e * A_c_s(t,MC))+ exp(e * A_nt_s(t,MC)));
omega_s_c(t,MC)=exp(e*A_c_s(t,MC))/(exp(e * A_f_s(t,MC)) + exp(e * A_c_s(t,MC))+ exp(e * A_nt_s(t,MC)));
omega_s_f(t,MC)=exp(e*A_nt_s(t,MC))/(exp(e * A_f_s(t,MC)) + exp(e * A_c_s(t,MC))+ exp(e * A_nt_s(t,MC)));
%% Expectations

Eq_tar_x(t,MC)=0;
Eq_sta_x(t,MC)=xq(t-1,MC);
Eq_ext_x(t,MC)=xq(t-1,MC)+alph*(xq(t-1,MC)-xq(t-64,MC));
Eq_tar_pi(t,MC)=0;
Eq_sta_pi(t,MC)=piq(t-1,MC);
Eq_ext_pi(t,MC)=piq(t-1,MC)+alph*(piq(t-1,MC)-piq(t-64,MC));
Et_s_c(t,MC)=s(t-1,MC)+kap_s*(s(t-1,MC)-s(t-2,MC));
Et_s_F(t,MC)=s(t-1,MC)+kap_s*(sf(t-1,MC)-s(t-2,MC));
    Eq_x(t,MC) =   omega_x_tar(t,MC) *Eq_tar_x(t,MC)+omega_x_sta(t,MC) *Eq_sta_x(t,MC)+omega_x_ext(t,MC) *Eq_ext_x(t,MC);  % market expectation of output gap
    Eq_pi(t,MC) =  omega_pi_tar(t,MC) *Eq_tar_pi(t,MC)+omega_pi_sta(t,MC) *Eq_sta_pi(t,MC)+omega_pi_ext(t,MC) *Eq_ext_pi(t,MC);   % market expectation of inflation
    E_s(t,MC)=omega_s_c(t,MC) *Et_s_c(t,MC)+omega_s_f(t,MC) *Et_s_F(t,MC); 
 
    %% demand of stock
    D_c(t,MC)=(Et_s_c(t,MC)-s(t-1,MC));
    D_f(t,MC)=(Et_s_F(t,MC)-s(t-1,MC));
    D(t,MC)=xq(t,MC)-c2*(s(t-1,MC)-pq(t,MC))-c3*i(t-1,MC);
%% State Space Representation
Y = [Eq_x(t,MC); Eq_pi(t,MC); E_s(t,MC)];
Z = [-h/sig *piq(t,MC)*delta_pi-h/sig *(delta_x-g*delta_s)*xq(t,MC);0; ...
    s(t-1,MC)*(1+a*k*c2+a*k*c3*delta_s)+a*omega_s_c(t,MC)*D_c(t,MC)+omega_s_f(t,MC)*D_f(t,MC)+eps_s(t,MC)];
Model = B*Y + Z;
X = A\Model;
x(t,MC)= X(1);
pi(t,MC)= X(2);
s(t,MC)=X(3);

%% A way to calculate yq and piq correctly is to preserve the daily
    % timeline, making sure that the value only changes once per quarter..
    if mod(t,64) == 0  %last day of a quarter 
        xq(t,MC) = mean(x(t-63:t,MC) );
        piq(t,MC) = mean( pi(t-63:t,MC) ); 
        pq(t,MC)=sum(piq(64:64:t,MC));

    else
        xq(t,MC) = xq(t-1,MC);
        piq(t,MC)= piq(t-1,MC);
        pq(t,MC)=pq(t-1,MC);

    end
    
    %% Taylor rule 
i(t,MC)= delta_pi*piq(t,MC)+delta_x*xq(t,MC)+delta_s*s(t,MC); %daily interest rate change only once quarterly if delta_s=0

% Generate quarterly results (with quarterly values ONLY)
xqq(1:T/64,MC) = xq(64:64:T,MC);
piqq(1:T/64,MC) = piq(64:64:T,MC);

%% Model 2

if t==129
    IMPULSE_Shock_var(t,MC)=IMPULSE_Shock;
else
    IMPULSE_Shock_var(t,MC)=0;

end
   eps_s2(t,MC)  = rho_s*eps_s2(t-1,MC) +sig_s*Randomz(t,MC) +IMPULSE_Shock_var(t,MC);

    sf2(t,MC)=g*xq2(t-1,MC);
 
 
 %% Attractivity
A_tar_pi2(t,MC) = zet*A_tar_pi2(t-1,MC) - (Eq_tar_pi2 (t-2,MC)-piq2(t-1,MC))^2; 
A_sta_pi2(t,MC) = zet*A_sta_pi2(t-1,MC) - (Eq_sta_pi2 (t-2,MC)-piq2(t-1,MC))^2; 
A_tar_x2(t,MC) = zet*A_tar_x2(t-1,MC) - (Eq_tar_x2 (t-2,MC)-xq2(t-1,MC))^2; 
A_sta_x2(t,MC) = zet*A_sta_x2(t-1,MC) - (Eq_sta_x2 (t-2,MC)-xq2(t-1,MC))^2; 
A_ext_x2(t,MC) = zet*A_ext_x2(t-1,MC) - (Eq_ext_x2 (t-2,MC)-xq2(t-1,MC))^2; 
A_ext_pi2(t,MC) = zet*A_ext_pi2(t-1,MC) - (Eq_ext_pi2 (t-2,MC)-piq2(t-1,MC))^2; 

A_f_s2(t,MC) = d*A_f_s2(t-1,MC) - (exp(s2 (t-2,MC))-exp(s2(t-1,MC)))*D_f2(t-2,MC); 
A_c_s2(t,MC) = d*A_c_s2(t-1,MC) - (exp(s2 (t-2,MC))-exp(s2(t-1,MC)))*D_c2(t-2,MC); 
A_nt_s2(t,MC)=0;

%% (((weights)))
omega_x_tar2(t,MC) = exp(phi*A_tar_x2(t,MC))/(exp(phi * A_tar_x2(t,MC)) + exp(phi * A_ext_x2(t,MC))+ exp(phi * A_sta_x2(t,MC)));
omega_pi_tar2(t,MC)  = exp(phi*A_tar_pi2(t,MC))/(exp(phi * A_tar_pi2(t,MC)) + exp(phi * A_ext_pi2(t,MC))+ exp(phi * A_sta_pi2(t,MC)));
omega_x_ext2(t,MC)  = exp(phi*A_ext_x2(t,MC))/(exp(phi * A_tar_x2(t,MC)) + exp(phi * A_ext_x2(t,MC))+ exp(phi * A_sta_x2(t,MC)));
omega_pi_ext2(t,MC)=exp(phi*A_ext_pi2(t,MC))/(exp(phi * A_tar_pi2(t,MC)) + exp(phi * A_ext_pi2(t,MC))+ exp(phi * A_sta_pi2(t,MC)));
omega_x_sta2(t,MC)=exp(phi*A_sta_x2(t,MC))/(exp(phi * A_tar_x2(t,MC)) + exp(phi * A_ext_x2(t,MC))+ exp(phi * A_sta_x2(t,MC)));
omega_pi_sta2(t,MC)=exp(phi*A_sta_pi2(t,MC))/(exp(phi * A_tar_pi2(t,MC)) + exp(phi * A_ext_pi2(t,MC))+ exp(phi * A_sta_pi2(t,MC)));

omega_s_nt2(t,MC)=exp(e*A_f_s2(t,MC))/(exp(e * A_f_s2(t,MC)) + exp(e * A_c_s2(t,MC))+ exp(e * A_nt_s2(t,MC)));
omega_s_c2(t,MC)=exp(e*A_c_s2(t,MC))/(exp(e * A_f_s2(t,MC)) + exp(e * A_c_s2(t,MC))+ exp(e * A_nt_s2(t,MC)));
omega_s_f2(t,MC)=exp(e*A_nt_s2(t,MC))/(exp(e * A_f_s2(t,MC)) + exp(e * A_c_s2(t,MC))+ exp(e * A_nt_s2(t,MC)));
%% Expectations

Eq_tar_x2(t,MC)=0;
Eq_sta_x2(t,MC)=xq2(t-1,MC);
Eq_ext_x2(t,MC)=xq2(t-1,MC)+alph*(xq2(t-1,MC)-xq2(t-64,MC));
Eq_tar_pi2(t,MC)=0;
Eq_sta_pi2(t,MC)=piq2(t-1,MC);
Eq_ext_pi2(t,MC)=piq2(t-1,MC)+alph*(piq2(t-1,MC)-piq2(t-64,MC));
Et_s_c2(t,MC)=s2(t-1,MC)+kap_s*(s2(t-1,MC)-s2(t-2,MC));
Et_s_F2(t,MC)=s2(t-1,MC)+kap_s*(sf2(t-1,MC)-s2(t-2,MC));
    Eq_x2(t,MC) =   omega_x_tar2(t,MC) *Eq_tar_x2(t,MC)+omega_x_sta2(t,MC) *Eq_sta_x2(t,MC)+omega_x_ext2(t,MC) *Eq_ext_x2(t,MC);  % market expectation of output gap
    Eq_pi2(t,MC) =  omega_pi_tar2(t,MC) *Eq_tar_pi2(t,MC)+omega_pi_sta2(t,MC) *Eq_sta_pi2(t,MC)+omega_pi_ext2(t,MC) *Eq_ext_pi2(t,MC);   % market expectation of inflation
    E_s2(t,MC)=omega_s_c2(t,MC) *Et_s_c2(t,MC)+omega_s_f2(t,MC) *Et_s_F2(t,MC); 
 
    %% demand of stock
    D_c2(t,MC)=(Et_s_c2(t,MC)-s2(t-1,MC));
    D_f2(t,MC)=(Et_s_F2(t,MC)-s2(t-1,MC));
    D2(t,MC)=xq2(t-64,MC)-c2*(s2(t-1,MC)-pq2(t-64,MC))-c3*i2(t-1,MC);
%% State Space Representation
Y2 = [Eq_x2(t,MC); Eq_pi2(t,MC); E_s2(t,MC)];
Z2 = [-h/sig *piq2(t,MC)*delta_pi-h/sig *(delta_x-g*delta_s)*xq2(t,MC);0; ...
    s2(t-1,MC)*(1+a*k*c2+a*k*c3*delta_s)+a*omega_s_c2(t,MC)*D_c2(t,MC)+omega_s_f2(t,MC)*D_f2(t,MC)+eps_s2(t,MC)];
Model2 = B*Y2 + Z2;
X2 = A\Model2;
x2(t,MC)= X2(1);
pi2(t,MC)= X2(2);
s2(t,MC)=X2(3);

%% A way to calculate yq and piq correctly is to preserve the daily
    % timeline, making sure that the value only changes once per quarter..
    if mod(t,64) == 0  %last day of a quarter 
        xq2(t,MC) = mean(x2(t-63:t,MC) );
        piq2(t,MC) = mean( pi2(t-63:t,MC) ); 
        pq2(t,MC)=sum(piq2(64:64:t,MC));

    else
        xq2(t,MC) = xq2(t-1,MC);
        piq2(t,MC)= piq2(t-1,MC);
        pq2(t,MC)=pq2(t-1,MC);


    end
    
    %% Taylor rule 
i2(t,MC)= delta_pi*piq2(t,MC)+delta_x*xq2(t,MC)+delta_s*s2(t,MC); %daily interest rate change only once quarterly if delta_s=0

% Generate quarterly results (with quarterly values ONLY)
xqq2(1:T/64,MC) = xq2(64:64:T,MC);
piqq2(1:T/64,MC) = piq2(64:64:T,MC);

    end
    end
%% Put all results into "result" variable
result.x  = mean(x2-x,2);
result.x1=x;
result.x2=x2;
result.x_irf = result.x (129:T);

result.xq = mean(xq2-xq,2);
result.xq_irf = result.xq (129:T);

result.pi= mean(pi2-pi,2);
result.pi_irf = result.pi (129:T);

result.piq = mean(piq2-piq,2);
result.piq_irf = result.piq (129:T);

result.i  = mean(i2-i,2);
result.i_irf = result.i (129:T);





result.Ex = mean(Et_x2-Et_x,2);
result.Ex_irf = result.Ex (129:T);

result.Es = mean(E_s2-E_s,2);
result.Es_irf = result.Es (129:T);

result.deltaEs = mean(E_s2-s2-E_s+s,2);
result.deltaEs_irf = result.deltaEs (129:T);



result.omega_x_tar = mean(omega_x_tar2-omega_x_tar,2);
result.omega_x_tar_irf = result.omega_x_tar (129:T);

result.omega_x_sta = mean(omega_x_sta2-omega_x_sta,2);
result.omega_x_sta_irf = result.omega_x_sta (129:T);

result.omega_x_ext = mean(omega_x_ext2-omega_x_ext,2);
result.omega_x_ext_irf = result.omega_x_ext (129:T);

result.omega_pi_ext = mean(omega_pi_ext2-omega_pi_ext,2);
result.omega_pi_ext_irf = result.omega_pi_ext (129:T);

result.Epi = mean(E_pi2-E_pi,2);
result.Epi_irf = result.Epi (129:T);

result.omega_pi_tar = mean(omega_pi_tar2-omega_pi_tar,2);
result.omega_pi_tar_irf = result.omega_pi_tar(129:T);

result.omega_pi_sta = mean(omega_pi_sta2-omega_pi_sta,2);
result.omega_pi_sta_irf = result.omega_pi_sta(129:T);

result.A_pi_tar=mean(A_tar_pi2-A_tar_pi,2);
result.A_pi_tar_irf = result.A_pi_tar(129:T);

result.A_pi_sta=mean(A_sta_pi2-A_sta_pi,2);
result.A_pi_sta_irf = result.A_pi_sta(129:T);

result.A_x_tar=mean(A_tar_x2-A_tar_x,2);
result.A_x_tar_irf = result.A_x_tar(129:T);

result.A_x_sta=mean(A_sta_x2-A_sta_x,2);
result.A_x_sta_irf = result.A_x_sta(129:T);

result.s=mean(s2-s,2);
result.s_irf = result.s (129:T);

result.IMPULSE_Shock_var = IMPULSE_Shock_var;
result.epss1 = eps_s;
result.epss2 = eps_s2;

result.omega_s_c=mean(omega_s_c2-omega_s_c,2);
result.omega_s_c_irf = result.omega_s_c (129:T);

result.omega_s_fund=mean(omega_s_f2-omega_s_f,2);
result.omega_s_fund_irf = result.omega_s_fund (129:T);

result.eps_s=mean(eps_s2-eps_s,2);
result.eps_s_irf = result.eps_s (129:T);

result.xqq = mean(xqq2-xqq,2);
result.xqq_irf = result.xqq (3:T/64);

result.piqq= mean(piqq2-piqq,2);
result.piqq_irf = result.piqq (3:T/64);






