function print_expectations(eqname, expectationmodelname, expectationmodelkind, withcalibration)

% Prints the exansion of the VAR_EXPECTATION or PAC_EXPECTATION term in files.
%
% INPUTS
% - eqname                      [string]    Name of the equation.
% - epxpectationmodelname       [string]    Name of the expectation model.
% - expectationmodelkind        [string]    Kind of the expectation model.
% - withcalibration             [logical]   Prints calibration if true.
%
% OUTPUTS
% None
%
% REMARKS
% print_expectations creates two text files
%
% - {expectationmodelname}-parameters.inc     which contains the declaration of the parameters specific to the expectation model kind term.
% - {expectationmodelname}-expression.inc     which contains the expanded version of the expectation model kind term.
%
% These routines are saved under the {modfilename}/model/{expectationmodelkind} subfolder, and can be
% used after in another mod file (ie included with the macro directive @#include).
%
% print_expectations also creates a matlab routine to evaluate the expectations (returning a dseries object).
%
% The variable expectationmodelkind can take two values 'var' or 'pac'.

% Copyright © 2018-2023 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare.  If not, see <https://www.gnu.org/licenses/>.

global M_

if nargin<4 || isempty(withcalibration)
    withcalibration = true;
end

% Check that the first input is a row character array.
if ~isrow(eqname)==1 || ~ischar(eqname)
    error('First input argument must be a row character array.')
end

% Check that the second input is a row character array.
if ~isrow(expectationmodelname)==1 || ~ischar(expectationmodelname)
    error('Second input argument must be a row character array.')
end

% Check that the third input is a row character array.
if ~isrow(expectationmodelkind)==1 || ~ischar(expectationmodelkind)
    error('Third input argument must be a row character array.')
end

% Check that the value of the second input is correct.
if ~ismember(expectationmodelkind, {'var', 'pac'})
    error('Wrong value for the second input argument.')
end

% Check that the model exists.
switch expectationmodelkind
  case 'var'
    if ~isfield(M_.var_expectation, expectationmodelname)
        error('VAR_EXPECTATION_MODEL %s is not defined.', expectationmodelname)
    else
        expectationmodelfield = 'var_expectation';
    end
  case 'pac'
    if ~isfield(M_.pac, expectationmodelname)
        error('PAC_EXPECTATION_MODEL %s is not defined.', expectationmodelname)
    else
        expectationmodelfield = 'pac';
    end
  otherwise
end

% Get the expectation model description
expectationmodel = M_.(expectationmodelfield).(expectationmodelname);

% Get the name of the associated VAR model and test its existence.
if isfield(expectationmodel, 'auxiliary_model_name') && ~isfield(M_.(expectationmodel.auxiliary_model_type), expectationmodel.auxiliary_model_name)
    switch expectationmodelkind
      case 'var'
        error('Unknown VAR/TREND_COMPONENT model (%s) in VAR_EXPECTATION_MODEL (%s)!', expectationmodel.auxiliary_model_name, expectationmodelname)
      case 'pac'
        error('Unknown VAR/TREND_COMPONENT model (%s) in PAC_EXPECTATION_MODEL (%s)!', expectationmodel.auxiliary_model_name, expectationmodelname)
      otherwise
    end
elseif isequal(expectationmodelkind, 'pac') && ~isfield(expectationmodel, 'auxiliary_model_name')
    error('print method does not work in PAC/MCE.')
end

auxmodel = M_.(expectationmodel.auxiliary_model_type).(expectationmodel.auxiliary_model_name);

%
% First print the list of parameters appearing in the VAR_EXPECTATION/PAC_EXPECTATION term.
%
if ~exist(sprintf('%s/model/%s', M_.fname, [expectationmodelkind '-expectations']), 'dir')
    mkdir(sprintf('%s/model/%s', M_.fname, [expectationmodelkind '-expectations']))
end

filename = sprintf('%s/model/%s/%s-parameters.inc', M_.fname, [expectationmodelkind '-expectations'], expectationmodelname);
fid = fopen(filename, 'w');
fprintf(fid, '// This file has been generated by dynare (%s).\n\n', datestr(now));

switch expectationmodelkind
  case 'var'
    parameter_declaration = 'parameters';
    for i=1:length(expectationmodel.param_indices)
        parameter_declaration = sprintf('%s %s', parameter_declaration, M_.param_names{expectationmodel.param_indices(i)});
    end
    fprintf(fid, '%s;\n\n', parameter_declaration);
    if withcalibration
        for i=1:length(expectationmodel.param_indices)
            fprintf(fid, '%s = %1.16f;\n', M_.param_names{expectationmodel.param_indices(i)}, M_.params(expectationmodel.param_indices(i)));
        end
    end
  case 'pac'
    parameter_declaration = 'parameters';
    if isfield(expectationmodel, 'h_param_indices')
        for i=1:length(expectationmodel.h_param_indices)
            parameter_declaration = sprintf('%s %s', parameter_declaration, M_.param_names{expectationmodel.h_param_indices(i)});
        end
    else
        for j=1:length(expectationmodel.components)
            for i=1:length(expectationmodel.components(j).h_param_indices)
                parameter_declaration = sprintf('%s %s', parameter_declaration, M_.param_names{expectationmodel.components(j).h_param_indices(i)});
            end
        end
    end
    fprintf(fid, '%s;\n\n', parameter_declaration);
    if withcalibration
        if isfield(expectationmodel, 'h_param_indices')
            for i=1:length(expectationmodel.h_param_indices)
                fprintf(fid, '%s = %1.16f;\n', M_.param_names{expectationmodel.h_param_indices(i)}, M_.params(expectationmodel.h_param_indices(i)));
            end
        else
            for j=1:length(expectationmodel.components)
                for i=1:length(expectationmodel.components(j).h_param_indices)
                    fprintf(fid, '%s = %1.16f;\n', M_.param_names{expectationmodel.components(j).h_param_indices(i)}, M_.params(expectationmodel.components(j).h_param_indices(i)));
                end
            end
        end
    end
    if isfield(expectationmodel, 'growth_neutrality_param_index')
        fprintf(fid, '\n');
        fprintf(fid, 'parameters %s;\n\n', M_.param_names{expectationmodel.growth_neutrality_param_index});
        if withcalibration
            fprintf(fid, '%s = %1.16f;\n', M_.param_names{expectationmodel.growth_neutrality_param_index}, M_.params(expectationmodel.growth_neutrality_param_index));
        end
        growth_correction = true;
    else
        growth_correction = false;
        if isfield(expectationmodel, 'components')
            for j=1:length(expectationmodel.components)
                if isfield(expectationmodel.components(j), 'growth_neutrality_param_index') && ~isempty(expectationmodel.components(j).growth_neutrality_param_index)
                    fprintf(fid, '\n');
                    fprintf(fid, 'parameters %s;\n\n', M_.param_names{expectationmodel.components(j).growth_neutrality_param_index});
                    if withcalibration
                        fprintf(fid, '%s = %1.16f;\n', M_.param_names{expectationmodel.components(j).growth_neutrality_param_index}, M_.params(expectationmodel.components(j).growth_neutrality_param_index));
                    end
                    growth_correction = true;
                end
            end
        end
    end
end
fclose(fid);

skipline()
fprintf('Parameters declarations and calibrations are saved in %s.\n', filename);

%
% Second print the expanded VAR_EXPECTATION/PAC_EXPECTATION term.
%

filename = sprintf('%s/model/%s/%s-expression.inc', M_.fname, [expectationmodelkind '-expectations'], expectationmodelname);
fid = fopen(filename, 'w');
fprintf(fid, '// This file has been generated by dynare (%s).\n', datestr(now));

switch expectationmodelkind
  case 'var'
    expression = write_expectations(expectationmodelname, expectationmodelkind, true);
  case 'pac'
    [expression, growthneutralitycorrection] = write_expectations(expectationmodelname, expectationmodelkind, true);
end

fprintf(fid, '%s', expression);
fclose(fid);

fprintf('Expectation unrolled expression is saved in %s.\n', filename);

%
% Second bis print the PAC growth neutrality correction term (if any).
%

if isequal(expectationmodelkind, 'pac') && growth_correction
    filename = sprintf('%s/model/%s/%s-growth-neutrality-correction.inc', M_.fname, [expectationmodelkind '-expectations'], expectationmodelname);
    fid = fopen(filename, 'w');
    fprintf(fid, '// This file has been generated by dynare (%s).\n', datestr(now));
    fprintf(fid, '%s', growthneutralitycorrection);
    fclose(fid);
    fprintf('Growth neutrality correction is saved in %s.\n', filename);
end

%
% Third print a routine for evaluating VAR_EXPECTATION/PAC_EXPECTATION term (returns a dseries object).
%
kind = [expectationmodelkind '_expectations'];
ndir = sprintf('+%s/+%s/+%s', M_.fname, kind, expectationmodelname);
if ~exist(ndir, 'dir')
    mkdir(sprintf('+%s/+%s/+%s', M_.fname, kind, expectationmodelname));
end
filename = sprintf('+%s/+%s/+%s/evaluate.m', M_.fname, kind, expectationmodelname);
fid = fopen(filename, 'w');

fprintf(fid, 'function ds = evaluate(dbase)\n\n');
fprintf(fid, '%% Evaluates %s term (%s).\n', kind, expectationmodelname);
fprintf(fid, '%%\n');
fprintf(fid, '%% INPUTS\n');
fprintf(fid, '%% - dbase     [dseries]  databse containing all the variables appearing in the auxiliary model for the expectation.\n');
fprintf(fid, '%%\n');
fprintf(fid, '%% OUTPUTS\n');
fprintf(fid, '%% - ds        [dseries]  the expectation term .\n');
fprintf(fid, '%%\n');
fprintf(fid, '%% REMARKS\n');
fprintf(fid, '%% The name of the appended variable in dbase is the declared name for the (PAC/VAR) expectation model.\n\n');
fprintf(fid, '%% This file has been generated by dynare (%s).\n\n', datestr(now));
fprintf(fid, 'ds = dseries();\n\n');

id = 0;

if isfield(expectationmodel, 'h_param_indices')
    decompose = false;
else
    if isequal(expectationmodelkind, 'pac')
        decompose = true;
    else
        decompose = false;
    end
end

clear('expression');

% Get coefficient values in the target (if any)
if exist(sprintf('+%s/pac_target_coefficients.m', M_.fname), 'file')
    targetcoefficients = feval(sprintf('%s.pac_target_coefficients', M_.fname), expectationmodelname, M_.params);
end

maxlag = max(auxmodel.max_lag);
if isequal(expectationmodel.auxiliary_model_type, 'trend_component')
    % Need to add a lag since the error correction equations are rewritten in levels.
    maxlag = maxlag+1;
end

if isequal(expectationmodelkind, 'var')
    timeindices = (0:(maxlag-1))+abs(expectationmodel.time_shift);
end

if isequal(expectationmodelkind, 'var') && isequal(expectationmodel.auxiliary_model_type, 'var')
    % Constant in the VAR auxiliary model
    id = id+1;
    expression = sprintf('%1.16f', M_.params(expectationmodel.param_indices(id)));
end

if isequal(expectationmodelkind, 'pac') && isequal(expectationmodel.auxiliary_model_type, 'var')
    % Constant in the VAR auxiliary model
    id = id+1;
    if isfield(expectationmodel, 'h_param_indices')
        constant = M_.params(expectationmodel.h_param_indices(id));
    else
        if decompose
            expressions = cell(length(expectationmodel.components), 1);
            for j=1:length(expectationmodel.components)
                expressions{j} = sprintf('%1.16f', M_.params(expectationmodel.components(j).h_param_indices(id)));
            end
        end
        constant = 0;
        for j=1:length(expectationmodel.components)
            constant = constant + targetcoefficients(j)*M_.params(expectationmodel.components(j).h_param_indices(id));
        end
    end
    if isfield(expectationmodel, 'h_param_indices')
        expression = sprintf('%1.16f', constant);
    end
end

for i=1:maxlag
    for j=1:length(auxmodel.list_of_variables_in_companion_var)
        id = id+1;
        variable = auxmodel.list_of_variables_in_companion_var{j};
        [variable, transformations] = rewrite_aux_variable(variable, M_);
        switch expectationmodelkind
          case 'var'
            parameter = M_.params(expectationmodel.param_indices(id));
          case 'pac'
            if isfield(expectationmodel, 'h_param_indices')
                parameter = M_.params(expectationmodel.h_param_indices(id));
            else
                parameter = 0;
                for k=1:length(expectationmodel.components)
                    parameter = parameter+targetcoefficients(k)*M_.params(expectationmodel.components(k).h_param_indices(id));
                end
            end
          otherwise
        end
        switch expectationmodelkind
          case 'var'
            if timeindices(i)>0
                variable = sprintf('dbase.%s(-%d)', variable, timeindices(i));
            else
                variable = sprintf('dbase.%s', variable);
            end
          case 'pac'
            variable = sprintf('dbase.%s(-%d)', variable, i);
          otherwise
        end
        if ~isempty(transformations)
            for k=length(transformations):-1:1
                variable = sprintf('%s.%s()', variable, transformations{k});
            end
        end
        if exist('expression','var')
            if parameter>=0
                expression = sprintf('%s+%1.16f*%s', expression, parameter, variable);
            elseif parameter<0
                expression = sprintf('%s-%1.16f*%s', expression, -parameter, variable);
            end
        else
            if parameter>=0
                expression = sprintf('%1.16f*%s', parameter, variable);
            elseif parameter<0
                expression = sprintf('-%1.16f*%s', -parameter, variable);
            end
        end
        if decompose
            for k=1:length(expectationmodel.components)
                parameter = M_.params(expectationmodel.components(k).h_param_indices(id));
                if parameter>=0
                    expressions{k} = sprintf('%s+%1.16f*%s', expressions{k}, parameter, variable);
                else
                    expressions{k} = sprintf('%s-%1.16f*%s', expressions{k}, -parameter, variable);
                end
            end
        end
    end
end

if isequal(expectationmodelkind, 'pac') && growth_correction
    if isfield(expectationmodel, 'growth_neutrality_param_index')
        pgrowth = M_.params(expectationmodel.growth_neutrality_param_index);
        for iter = 1:numel(expectationmodel.growth_linear_comb)
            vgrowth='';
            variable = [];
            if expectationmodel.growth_linear_comb(iter).exo_id > 0
                variable = M_.exo_names{expectationmodel.growth_linear_comb(iter).exo_id};
            elseif expectationmodel.growth_linear_comb(iter).endo_id > 0
                variable = M_.endo_names{expectationmodel.growth_linear_comb(iter).endo_id};
            end
            if ~isempty(variable)
                [variable, transformations] = rewrite_aux_variable(variable, M_);
                if isempty(transformations)
                    if expectationmodel.growth_linear_comb(iter).lag ~= 0
                        variable = sprintf('%s(%d)', variable, expectationmodel.growth_linear_comb(iter).lag);
                    end
                else
                    for k=rows(transformations):-1:1
                        if isequal(transformations{k,1}, 'lag')
                            variable = sprintf('%s.lag(%u)', variable, -transformations{k,2});
                        elseif isequal(transformations{k,1}, 'diff')
                            if isempty(transformations{k,2})
                                variable = sprintf('%s.diff()', variable);
                            else
                                variable = sprintf('%s.lag(%u).diff()', variable, transformations{k,2});
                            end
                        else
                            variable = sprintf('%s.%s()', variable, transformations{k});
                        end
                    end
                end
                vgrowth = strcat('dbase.', variable);
            end
            if expectationmodel.growth_linear_comb(iter).param_id > 0
                if ~isempty(vgrowth)
                    vgrowth = sprintf('%1.16f*%s',M_.params(expectationmodel.growth_linear_comb(iter).param_id), vgrowth);
                else
                    vgrowth = num2str(M_.params(expectationmodel.growth_linear_comb(iter).param_id), '%1.16f');
                end
            end
            if abs(expectationmodel.growth_linear_comb(iter).constant) ~= 1
                if ~isempty(vgrowth)
                    vgrowth = sprintf('%1.16f*%s', expectationmodel.growth_linear_comb(iter).constant, vgrowth);
                else
                    vgrowth = num2str(expectationmodel.growth_linear_comb(iter).constant, '%1.16f');
                end
            end
            if iter > 1
                if expectationmodel.growth_linear_comb(iter).constant > 0
                    linearCombination = sprintf('%s+%s', linearCombination, vgrowth);
                else
                    linearCombination = sprintf('%s-%s', linearCombination, vgrowth);
                end
            else
                    linearCombination = vgrowth;
                end
        end % loop over growth linear combination elements
        growthcorrection = sprintf('%1.16f*(%s)', pgrowth, linearCombination);
    else
        first = true;
        for i=1:length(expectationmodel.components)
            if ~isequal(expectationmodel.components(i).kind, 'll') && isfield(expectationmodel.components(i), 'growth_neutrality_param_index') && isfield(expectationmodel.components(i), 'growth_linear_comb') && ~isempty(expectationmodel.components(i).growth_linear_comb)
                pgrowth = targetcoefficients(i)*M_.params(expectationmodel.components(i).growth_neutrality_param_index);
                for iter = 1:numel(expectationmodel.components(i).growth_linear_comb)
                    vgrowth='';
                    variable=[];
                    if expectationmodel.components(i).growth_linear_comb(iter).exo_id > 0
                        variable = M_.exo_names{expectationmodel.components(i).growth_linear_comb(iter).exo_id};
                    elseif expectationmodel.components(i).growth_linear_comb(iter).endo_id > 0
                        variable = M_.endo_names{expectationmodel.components(i).growth_linear_comb(iter).endo_id};
                    end
                    if ~isempty(variable)
                        [variable, transformations] = rewrite_aux_variable(variable, M_);
                        if isempty(transformations)
                            if expectationmodel.components(i).growth_linear_comb(iter).lag ~= 0
                                variable = sprintf('%s(%d)', variable, expectationmodel.components(i).growth_linear_comb(iter).lag);
                            end
                        else
                            for k=rows(transformations):-1:1
                                if isequal(transformations{k,1}, 'lag')
                                    variable = sprintf('%s.lag(%u)', variable, -transformations{k,2});
                                elseif isequal(transformations{k,1}, 'diff')
                                    if isempty(transformations{k,2})
                                        variable = sprintf('%s.diff()', variable);
                                    else
                                        variable = sprintf('%s.lag(%u).diff()', variable, transformations{k,2});
                                    end
                                else
                                    variable = sprintf('%s.%s()', variable, transformations{k});
                                end
                            end
                        end
                        vgrowth = strcat('dbase.', variable);
                    end
                    if expectationmodel.components(i).growth_linear_comb(iter).param_id > 0
                        if ~isempty(vgrowth)
                            vgrowth = sprintf('%1.16f*%s',M_.params(expectationmodel.components(i).growth_linear_comb(iter).param_id), vgrowth);
                        else
                            vgrowth = num2str(M_.params(expectationmodel.components(i).growth_linear_comb(iter).param_id), '%1.16f');
                        end
                    end
                    if abs(expectationmodel.components(i).growth_linear_comb(iter).constant) ~= 1
                        if ~isempty(vgrowth)
                            vgrowth = sprintf('%1.16f*%s', expectationmodel.components(i).growth_linear_comb(iter).constant, vgrowth);
                        else
                            vgrowth = num2str(expectationmodel.components(i).growth_linear_comb(iter).constant, '%1.16f');
                        end
                    end
                    if iter > 1
                        if expectationmodel.components(i).growth_linear_comb(iter).constant > 0
                            linearCombination = sprintf('%s+%s', linearCombination, vgrowth);
                        else
                            linearCombination = sprintf('%s-%s', linearCombination, vgrowth);
                        end
                    else
                        linearCombination = vgrowth;
                    end
                end % loop over growth linear combination elements
                if first
                    growthcorrection = sprintf('%1.16f*(%s)', pgrowth, linearCombination);
                    first = false;
                else
                    if pgrowth>0
                        growthcorrection = sprintf('%s+%1.16f*(%s)', growthcorrection, pgrowth, linearCombination);
                    elseif pgrowth<0
                        growthcorrection = sprintf('%s-%1.16f*(%s)', growthcorrection, -pgrowth, linearCombination);
                    end
                end
            end
        end
    end
    expression = sprintf('%s+%s', expression, growthcorrection);
end % growth_correction

fprintf(fid, 'ds.%s = %s;\n', expectationmodelname, expression);
if exist('expressions', 'var')
    for i=1:length(expressions)
        fprintf(fid, 'ds.%s = %s;\n', M_.lhs{expectationmodel.components(i).aux_id}, expressions{i});
    end
end
fclose(fid);

fprintf('Expectation dseries expression is saved in %s.\n', filename);

skipline();

rehash