/* This template is used to produce IRFs
**
** Requires the following definitions:
**
** - @#define ShockName             [string]     Name of the exogenous variables.
** - @#define ShockValue            [float]      scalar, value of the exogenous variable in first period.
** - @#define VariableName          [string]     Name of the endogenous variables.
** - @#define VariablePdev          [string]     Endogenous variables in dev (see below).
** - @#define VariablePpdev         [string]     Endogenous variables in pdev (see below).
** - @#define VariablePqdev         [string]     Endogenous variables in qdev (see below).
** - @#define VariablePadev         [string]     Endogenous variables in adev (see below).
** - @#define VariablePydev         [string]     Endogenous variables in ydev (see below).
** - @#define SimulationPeriods     [integer]    scalar, number of periods for the perfect foresight solver.
** - @#define IrfPeriods            [integer]    scalar, number of periods for the IRFs
**
** 2020-10-13 stepan@adjemian.eu
*/

@#include "update-database.mod"

// Check that the baseline is fine (solver should exit on first iteration)
initval_file(series=baseline, first_obs=3);
perfect_foresight_setup(periods=@{SimulationPeriods});
perfect_foresight_solver(stack_solve_algo=0);

Baseline = dseries([transpose(oo_.endo_simul), oo_.exo_simul], 5000Q1, [M_.endo_names; M_.exo_names]);

variables = [dev, pdev, adev, qdev, ydev];

// Subplots parameters
nvar = length(variables); // Number of subplots is the number of endogenous variables
nrow = ceil(sqrt(nvar));  // Number of rows
ncol = ceil(nvar/nrow);   // Number of columns

disp('*****************************************************')
disp('               @{ShockName}                          ')
disp('*****************************************************')

shocks;
  var @{ShockName};
  periods 1;
  values @{ShockValue};
end;

initval_file(series=baseline, first_obs=3);
perfect_foresight_setup(periods=@{SimulationPeriods});
perfect_foresight_solver(stack_solve_algo=0, no_homotopy);

irfs = dseries('./mat/irfs_backward_model_@{ShockName}.mat');

Simulated_time_series = dseries(transpose(oo_.endo_simul), Baseline.dates(1), M_.endo_names);
Simulated_time_series.save(sprintf('./mat/%s_orig_irf_%s', M_.fname, '@{ShockName}'), 'mat');
Simulated_time_series = irf_trans_report(Baseline(Baseline.dates(1):Baseline.dates(end)), Simulated_time_series, 'dev', dev, 'pdev', pdev, 'qdev', qdev, 'adev', adev, 'ydev', ydev, 'start_date', Baseline.dates(1));
irfs = irf_trans_report(baseline, irfs, 'dev', dev, 'pdev', pdev, 'qdev', qdev, 'adev', adev, 'ydev', ydev, 'start_date', irfs.dates(1));

@#if BUILD_MATLAB_PLOTS

    if ~isfolder('./pdf')
        mkdir('./pdf')
    end

    figure('Name', sprintf('Shock to @{ShockName} (@{ShockValue})'), 'visible', 'off');
    for i=1:nvar
        subplot(nrow, ncol, i)
        plot(0:120, irfs{variables{i}}.data(1:@{IrfPeriods}+1), '-k', 'linewidth', .8)
        hold on
        plot(0:120, Simulated_time_series{variables{i}}.data(1:@{IrfPeriods}+1), '-', 'linewidth', .8)
        hold off
        axis tight
        title(variables{i}, 'Interpreter', 'none')
    end

    set(findall(gcf,'-property','FontSize'),'FontSize',5)
    set(findall(gcf, 'Type', 'Text'),'FontWeight', 'Normal')
    print(sprintf('./pdf/%s_irf_%s.pdf', M_.fname, '@{ShockName}'), '-dpdf', '-bestfit');

@#endif

Simulated_time_series.save(sprintf('./mat/%s_irf_%s', M_.fname, '@{ShockName}'), 'mat');

@#ifdef SUBFOLDERNAME

    % Save data for PGF plots
    for i=1:length(variables)
        if ~exist('../irfs-data/backward/@{ShockName}')
            mkdir('../irfs-data/backward/@{ShockName}');
        end
        if ~exist(sprintf('../irfs-data/hybrid/%s/@{ShockName}', '@{SUBFOLDERNAME}'))
            mkdir(sprintf('../irfs-data/hybrid/%s/@{ShockName}', '@{SUBFOLDERNAME}'));
        end
        fid = fopen(sprintf('../irfs-data/backward/@{ShockName}/%s.dat', variables{i}), 'w');
        fprintf(fid, '%d \t %12.6f\n', [0:@{IrfPeriods}; transpose(irfs{variables{i}}.data(1:@{IrfPeriods}+1))]);
        fclose(fid);
        fid = fopen(sprintf('../irfs-data/hybrid/%s/@{ShockName}/%s.dat', '@{SUBFOLDERNAME}', variables{i}), 'w');
        fprintf(fid, '%d \t %12.6f\n', [0:@{IrfPeriods}; transpose(Simulated_time_series{variables{i}}.data(1:@{IrfPeriods}+1))]);
        fclose(fid);
    end

@#endif

shocks(overwrite);
end;

initval_file(series=baseline, first_obs=3);
perfect_foresight_setup(periods=@{SimulationPeriods});
perfect_foresight_solver(stack_solve_algo=0);
