1function innovation_paths = reversed_extended_path(controlled_variable_names, control_innovation_names, dataset)
2% Inversion of the extended path simulation approach. This routine computes the innovations needed to
3% reproduce the time path of a subset of endogenous variables. The initial condition is teh deterministic
4% steady state.
5%
6% INPUTS
7%  o controlled_variable_names        [string]    n*1 matlab's cell.
8%  o control_innovation_names         [string]    n*1 matlab's cell.
9%  o dataset                          [structure]
10% OUTPUTS
11%  o innovations                      [double]  n*T matrix.
12%
13% ALGORITHM
14%
15% SPECIAL REQUIREMENTS
16
17% Copyright (C) 2010-2018 Dynare Team.
18%
19% This file is part of Dynare.
20%
21% Dynare is free software: you can redistribute it and/or modify
22% it under the terms of the GNU General Public License as published by
23% the Free Software Foundation, either version 3 of the License, or
24% (at your option) any later version.
25%
26% Dynare is distributed in the hope that it will be useful,
27% but WITHOUT ANY WARRANTY; without even the implied warranty of
28% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
29% GNU General Public License for more details.
30%
31% You should have received a copy of the GNU General Public License
32% along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
33
34global M_ oo_ options_
35
36%% Initialization
37
38% Load data.
39eval(dataset.name);
40dataset.data = [];
41for v = 1:dataset.number_of_observed_variables
42    eval(['dataset.data = [ dataset.data , ' dataset.variables(v,:) ' ];'])
43end
44data = dataset.data(dataset.first_observation:dataset.first_observation+dataset.number_of_observations,:);
45
46% Compute the deterministic steady state.
47steady_;
48
49%  Compute the first order perturbation reduced form.
50old_options_order = options_.order; options_.order = 1;
51[dr,info,M_,options_,oo_] = resol(0,M_,options_,oo_);
52oo_.dr = dr;
53options_.order = old_options_order;
54
55% Set various options.
56options_.periods = 100;
57
58% Set-up oo_.exo_simul.
59oo_=make_ex_(M_,options_,oo_);
60
61% Set-up oo_.endo_simul.
62oo_=make_y_(M_,options_,oo_);
63
64% Get indices of the controlled endogenous variables in endo_simul.
65n  = length(controlled_variable_names);
66iy = NaN(n,1);
67for k=1:n
68    iy(k) = strmatch(controlled_variable_names{k}, M_.endo_names, 'exact');
69end
70
71% Get indices of the controlled endogenous variables in dataset.
72iy_ = NaN(n,1);
73for k=1:n
74    iy_(k) = strmatch(controlled_variable_names{k},dataset.variables,'exact');
75end
76
77% Get indices of the control innovations in exo_simul.
78ix = NaN(n,1);
79for k=1:n
80    ix(k) = strmatch(control_innovation_names{k},M_.exo_names,'exact');
81end
82
83% Get the length of the sample.
84T = size(data,1);
85
86% Output initialization.
87innovation_paths = zeros(n,T);
88
89% Initialization of the perfect foresight model solver.
90perfect_foresight_simulation();
91
92% Set options for fsolve.
93options = optimset('MaxIter',10000,'Display','Iter');
94
95%% Call fsolve recursively
96for t=1:T
97    x0 = zeros(n,1);
98    y_target  = transpose(data(t,iy_));
99    total_variation = y_target-transpose(oo_.endo_simul(t+M_.maximum_lag,iy));
100    for i=1:100
101        [t,i]
102        y = transpose(oo_.endo_simul(t+M_.maximum_lag,iy)) + (i/100)*y_target
103        [tmp,fval,exitflag] = fsolve('ep_residuals', x0, options, y, ix, iy, oo_.steady_state, oo_.dr, M_.maximum_lag, M_.endo_nbr);
104    end
105    if exitflag==1
106        innovation_paths(:,t) = tmp;
107    end
108    % Update endo_simul.
109    oo_.endo_simul(:,1:end-1) = oo_.endo_simul(:,2:end);
110    oo_.endo_simul(:,end) = oo_.steady_state;
111end