// --+ options: json=compute +--

@#include "core.mod"

@#if !BACKWARD

    dev = {};
    @#ifdef VariablePdev
        @#if length(VariablePdev)>0
            @#for i in 1:length(VariablePdev)
                dev{@{i}} = '@{VariablePdev[i]}';
            @#endfor
        @#endif
    @#endif

    pdev = {};
    @#ifdef VariablePpdev
        @#if length(VariablePpdev)>0
            @#for i in 1:length(VariablePpdev)
                pdev{@{i}} = '@{VariablePpdev[i]}';
            @#endfor
        @#endif
    @#endif

    qdev = {};
    @#ifdef VariablePqdev
        @#if length(VariablePqdev)>0
            @#for i in 1:length(VariablePqdev)
                qdev{@{i}} = '@{VariablePqdev[i]}';
            @#endfor
        @#endif
    @#endif

    adev = {};
    @#ifdef VariablePadev
        @#if length(VariablePadev)>0
            @#for i in 1:length(VariablePadev)
                adev{@{i}} = '@{VariablePadev[i]}';
            @#endfor
        @#endif
    @#endif

    ydev = {};
    @#ifdef VariablePydev
        @#if length(VariablePydev)>0
            @#for i in 1:length(VariablePydev)
                ydev{@{i}} = '@{VariablePydev[i]}';
            @#endfor
        @#endif
    @#endif

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

    @#if BUILD_IRFS

        @#define SimulationPeriods  = 1500

        @#ifndef IrfPeriods
        @#define IrfPeriods = 120
        @#endif

        @#if defined(ShockValues) && !defined(ShockNames)
            @#error "Definition of ShockNames is missing"
        @#endif

        @#ifdef ShockNames
            @#if !isarray(ShockNames)
                @#error "ShockNames must be an array"
            @#endif
            @#ifndef ShockValues
                @#error "Definition of ShockValues is missing"
            @#endif
            @#if !isarray(ShockValues)
                @#error "ShockValues must be an array"
            @#endif
            @#if !(length(ShockNames)==length(ShockValues))
                @#error "ShockValues and ShockNames arrays must have the same number of elements"
            @#endif
            @#for i in 1:length(ShockNames)
                @#define ShockName = ShockNames[i]
                @#define ShockValue = ShockValues[i]
                @#if !isstring(ShockName)
                    @#error "ShockNames[@{i}] has to be a string (name of an exogenous variable)"
                @#endif
                @#if !isreal(ShockValue)
                     @#error "ShockValues[@{i}] has to be a real number"
                @#endif
                @#include "./templates/build-irfs.mod"
            @#endfor
        @#else
            @#ifndef ShockName
                @#error "Define ShockName and ShockValue or ShockNames and ShockValues"
            @#endif
            @#ifndef ShockValue
                @#error "Define ShockName and ShockValue or ShockNames and ShockValues"
            @#endif
            @#if !isstring(ShockName)
                @#error "ShockName has to be a string (name of an exogenous variable)"
            @#endif
            @#if !isreal(ShockValue)
                @#error "ShockValue has to be a real number"
            @#endif
            @#include "./templates/build-irfs.mod"
        @#endif

    @#elseif FORWARD_GUIDANCE_PUZZLE

        @#include "update-database.mod"

        initval_file(series=baseline, first_obs=3);
        perfect_foresight_setup(periods=1500);
        perfect_foresight_solver(stack_solve_algo=0);

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

        nperiods = options_.periods;

        irfs = cell(20,1);

        @#for i in [i for i in 1:20]

            shocks;
                var res_U2_STN;
                periods @{i};
                values 1;
            end;

            initval_file(series=baseline, first_obs=3);
            perfect_foresight_setup(periods=1500);
            perfect_foresight_solver(stack_solve_algo=0, no_homotopy);

            irfs{@{i}} = dseries(Simulated_time_series.data, Baseline.dates(5), Simulated_time_series.name);

            shocks(overwrite);
            end;

            initval_file(series=baseline, first_obs=3);
            perfect_foresight_setup(periods=1500);
            perfect_foresight_solver(stack_solve_algo=0);

        @#endfor

        save(sprintf('irfs_data_fgp_%s', M_.fname), 'irfs');

        @#if BUILD_MATLAB_PLOTS

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

            ColourMatrix = [0 0 1;
                1 0 0];

            figure;

            % U2_STN
            s1 = subplot(2, 2, 1);
            hold on
            for i=1:20
                lambda = 1-i/20;
                plot(irfs{i}.U2_STN.data(1:60), 'linewidth', 1.5, 'color', lambda*ColourMatrix(1,:) + (1-lambda)*ColourMatrix(2,:))
            end
            hold off
            grid on
            box on
            title('U2_STN','Interpreter','none')
            set(s1,'FontName','cmr12');

            % U2_G_YER
            s2 = subplot(2, 2, 2);
            hold on
            for i=1:20
                lambda = 1-i/20;
                plot(irfs{i}.U2_G_YER.data(1:60), 'linewidth', 1.5, 'color', lambda*ColourMatrix(1,:) + (1-lambda)*ColourMatrix(2,:))
            end
            hold off
            grid on
            box on
            title('U2_G_YER','Interpreter','none')
            set(s2,'FontName','cmr12');

            % U2_A_YED
            s3 = subplot(2, 2, 3);
            hold on
            for i=1:20
                lambda = 1-i/20;
                plot(irfs{i}.U2_A_YED.data(1:60), 'linewidth', 1.5, 'color', lambda*ColourMatrix(1,:) + (1-lambda)*ColourMatrix(2,:))
            end
            hold off
            grid on
            box on
            title('U2_A_YED','Interpreter','none')
            set(s3,'FontName','cmr12');

        @#endif

        @#ifdef SUBFOLDERNAME
            % Save data for PGF plots
            for i=1:length(variables)
                if ~exist(sprintf('../irfs-data/hybrid/%s', '@{SUBFOLDERNAME}'))
                    mkdir(sprintf('../irfs-data/hybrid/%s', '@{SUBFOLDERNAME}'));
                end
                fid = fopen(sprintf('../irfs-data/hybrid/%s/%s.dat', '@{SUBFOLDERNAME}', variables{i}), 'w');
                fprintf(fid, '%d \t %12.6f \t %12.6f \t %12.6f \t %12.6f \t %12.6f \t %12.6f \t %12.6f \t %12.6f \t %12.6f \t %12.6f \t %12.6f \t %12.6f \t %12.6f \t %12.6f \t %12.6f \t %12.6f \t %12.6f \t %12.6f \t %12.6f \t %12.6f\n', [1:60; transpose(irfs{1}.(variables{i}).data(1:60)); transpose(irfs{2}.(variables{i}).data(1:60)); transpose(irfs{3}.(variables{i}).data(1:60)); transpose(irfs{4}.(variables{i}).data(1:60)); transpose(irfs{5}.(variables{i}).data(1:60)); transpose(irfs{6}.(variables{i}).data(1:60)); transpose(irfs{7}.(variables{i}).data(1:60)); transpose(irfs{8}.(variables{i}).data(1:60)); transpose(irfs{9}.(variables{i}).data(1:60)); transpose(irfs{10}.(variables{i}).data(1:60)); transpose(irfs{11}.(variables{i}).data(1:60)); transpose(irfs{12}.(variables{i}).data(1:60)); transpose(irfs{13}.(variables{i}).data(1:60)); transpose(irfs{14}.(variables{i}).data(1:60)); transpose(irfs{15}.(variables{i}).data(1:60)); transpose(irfs{16}.(variables{i}).data(1:60)); transpose(irfs{17}.(variables{i}).data(1:60)); transpose(irfs{18}.(variables{i}).data(1:60)); transpose(irfs{19}.(variables{i}).data(1:60)); transpose(irfs{20}.(variables{i}).data(1:60))]);
                fclose(fid);
            end
       @#endif

       @#if BUILD_MATLAB_PLOTS
           lgd = legend('t=1', 't=2', 't=3', 't=4', 't=5', 't=6', 't=7', 't=8', 't=9', 't=10', ...
                         't=11', 't=12', 't=13', 't=14', 't=15', 't=16', 't=17', 't=18', 't=19', 't=20');

            % Programatically move the Legend
            set(lgd, ...
                'Position', [0.634934617934261 0.0803022281216621 0.202486960333143 0.380774024125193], ...
                'NumColumns', 2, ...
                'LineWidth', 1, ...
                'FontSize', 18, ...
                'FontName', 'cmr12', ...
                'Box', 'off');

            print(sprintf('./eps/fgp_%s', M_.fname), '-depsc2');

       @#endif

    @#elseif EXPECTED_PRODUCTIVITY_SHOCK

        @#include "update-database.mod"

        initval_file(series=baseline, first_obs=3);
        perfect_foresight_setup(periods=1500);
        perfect_foresight_solver(stack_solve_algo=0);

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

        nperiods = options_.periods;

        // 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

        // Set the number of periods.
        speriods = 120;

        shockvalues = [0.01 0.01 0.01 0.01];// 0.01 0.01 0.01 0.01];

        shocks;
            var res_U2_Q_T_PRO;
            periods 44:47;
            values (shockvalues);
        end;

        initval_file(series=baseline, first_obs=3);
        perfect_foresight_setup(periods=1500);
        perfect_foresight_solver(stack_solve_algo=0, no_homotopy);

        irfs = dseries('mat/irfs_backward_model_res_U2_Q_T_PRO.mat');
        Simulated_time_series = dseries(Simulated_time_series.data, Baseline.dates(1), Simulated_time_series.name);
        Simulated_time_series = irf_trans_report(Baseline, 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', 'Shock to res_U2_Q_T_PRO (.01)', 'visible', 'off');
            for i=1:nvar
                subplot(nrow, ncol, i)
                plot(Simulated_time_series{variables{i}}.data(2:(2+speriods)), '-', 'linewidth', 2)
                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/irf_%s.pdf', 'res_U2_Q_T_PRO'),'-dpdf', '-bestfit');

        @#endif

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

    @#endif

@#endif

if isfile('./baseline_updated.mat')
    delete baseline_updated.mat
end

if isfile('./baseline_updated.old.mat')
    delete baseline_updated.old.mat
end

if isfile('./hybrid.log')
    delete hybrid.log
end

if isfile('./hybrid_results.mat')
    delete hybrid_results.mat
end

if isfolder('./+hybrid')
    rmdir('+hybrid', 's');
    rmdir('hybrid', 's');
end
