%*************************************************************************
% MAIN_PROG - Generate the data (including model impulse-responses)
% for "A Model-Based Evaluation of the Debate on the Size of the Tax 
% Multiplier" (Chahrour, Schmitt-Grohé, Uribe; American Economic Journal: 
% Economic Policy; 2011)
%*************************************************************************


%Simulation parameters
nk    = 1000;          %# of replications
ndrop = 750;           %# burn in periods in each simulation
nobs  = 250;          %# of observations in each simulation
nper  = ndrop + nobs;  %Total # of periods simulated

%Parameters for output
nper_ir = 20;    %# of periods in impulse responses
ny_ir   = 3;     %# of observables
nlag_ir = 4;     %# of lags in BP specification
nl      = 12;    %# of lags in RR specification 
neps    = 8;     %# of fundamental shocks
 
%Include constant in regressions?
incl_c = 0;   

%Initialize some matrices
ir_bp  = zeros(nper_ir,ny_ir,nk,2);
ir_fg0 = zeros(nper_ir,ny_ir,nk,2);
ir_rr  = zeros(nper_ir,ny_ir,nk,2);
ir_rrg = zeros(nper_ir,ny_ir,nk,2);

%Main Loop: 1) Solve Model  2) Simulate Data  3) Compute IR's
for ll = 1:2
  
    %Get derivatives of model FOC's
    if ll == 1
        %Baseline Case
        [fy fx fyp fxp eta param] = model(parameters_baseline);
    elseif ll == 2
        %Anticipation Case
        [fy fx fyp fxp eta param] = model(parameters_anticipation);
    end

    %Compute model transition and observation matrices
    [gx,hx]=gx_hx(fy,fx,fyp,fxp);

    %Save model solution for figure 5
    save('figures_and_tables/model_solution', 'gx', 'hx', 'eta');
    
    %Parameters used for scaling responses
    tbar = param.tbar;
    ybar = param.ybar;
    sig  = param.sig;
    
    %Unit shock to model
    mu = zeros(neps,1);
    mu(1) = 1;
    
    %Normalize initial shock to give unit theory tax response
    etat = eta;
    etat(eta~=0) = 1;
    scl = gx*etat*mu*tbar/ybar;
    scl = scl(15);
    x0 = etat*mu/scl;
    
    %Compute model impulse responses 
    [waste irY irX] = ir(gx,hx, x0,nper_ir);  %Normalizing IR
    ir_mod = [irY(:,15)*tbar/ybar, irY(:,16)*sig, irY(:,9)];
    
    %Impose identification assumption for BP:
    %shocks are correctly identified => impact response is correct
    ir0 = ir_mod(1,:);

    %Scale of initial shock for RR:
    z0 = 1/scl;     
    
    %Simulate and compute multipliers
    for k = 1:nk
        
        %Simulate Data
        [sx, sy, seps] = sim_dat(gx,hx,eta,nper,ndrop);
        
        %Relevant Observables: [TAX/Y,G,Y]
        Y = [sy(15,:)'*tbar/ybar, sy(16,:)'*sig, sy(9,:)'];
        
        %Realized Tax Shocks
        eps_tax = sx(14,:)';           %Surprise Tax Shock
        eps_tax_news = sx(15,:)';      %News shock implemented today
        news_shocks = sx(16:16+5,:)';  %News about taxes in the next six periods
        
        %Run Blanchard-Perotti (2002) VAR
        [bet] = quick_var(Y, nlag_ir,incl_c);
        ir_bp(:,:,k,ll) = quick_ir(bet, nper_ir, ir0);
        
        %Run Favero-Giavazzi (2010), with controls that are only 1 period long
        [bet, bet1, bet2] = quick_var(Y, nlag_ir,incl_c,eps_tax+eps_tax_news,0);
        ir_fg0(:,:,k,ll)  = quick_ir(bet1, nper_ir, zeros(1,3) ,bet2, z0);
     
        for j = 1:ny_ir
                 
            %Romer-Romer (2007) in levels
            [bet, bet1, bet2] = quick_var(Y(:,j),0,incl_c,eps_tax+eps_tax_news,nl); 
            ir_rr(:,j,k,ll)   = quick_ir(bet1, nper_ir, 0, bet2, z0);
            
            %Romer-Romer in log-changes
            [bet, bet1, bet2] = quick_var(diff(Y(:,j)),0,incl_c,eps_tax(2:end)+eps_tax_news(2:end),nl);
            ir_rrg(:,j,k,ll)  = cumsum(quick_ir(bet1, nper_ir, 0, bet2, z0));
           
        end
      
        %Display Counter
        if mod(k,100)==0
            disp(num2str(k));
        end
    end
end


%Save required data
mean_bp         = mean(ir_bp(:,:,:,1),3);
mean_bp(:,:,2)  = mean(ir_bp(:,:,:,2),3);

mean_rr         = mean(ir_rr(:,:,:,1),3);
mean_rr(:,:,2)  = mean(ir_rr(:,:,:,2),3);

mean_rrg        = mean(ir_rrg(:,:,:,1),3);
mean_rrg(:,:,2) = mean(ir_rrg(:,:,:,2),3);

mean_fg0        = mean(ir_fg0(:,:,:,1),3);
mean_fg0(:,:,2) = mean(ir_fg0(:,:,:,2),3); 

%save(['figures_and_tables/figures_data_' num2str(nobs)], 'ir_bp', 'ir_mod', 'ir_rr', 'mean_bp', 'mean_fg0', 'mean_rr', 'mean_rrg')

% Uncomment to make figures automatically
% cd figures_and_tables/
% make_figures
% cd ..
