function [M_,oo_]=stoch_simul_para_t(opt, M_orig,options_orig,oo_orig)

% This function computes the solution for a linearized model with the
% anticipated parameter changes, as explained in the companion document.
% The input is a structure with the following fields:
%
% names_change = names of parameters that change, as declared in dynare.
%              The parameter SOLVE should NOT be included. 
%              All these need to have been calibrated with the 
%              'set_param_value' option in the mod-file.
%              Define the number of these parameter as N_change.
%           T0 = moment of the announcement, >=1.
%        Tstar = moment of the terminal point >T0.
%         Tmax = Number of periods to be considered (with Tmax>=Tstar)
% values_theta = values for the parameters that are assumed to change. If
%              the vector of parameter is assumed to be constant between
%              the periods T0 and Tstar, this should be a matrix of size
%              N_change by 3 (each row ordered as the declaration in 
%              names_change). The first column shoulddyrjs contain the values for
%              t<T0. The second for T0<=t<Tstar. The third for t>=Tstar.
%              If a given sequence for the parameters is desired, this
%              should be a matrix of size N_change by 3N_change by
%              (2+Tstar-T0). The first column should contain the values for
%              t<T0. from column 2 to Tstar-T0+1, the sequence for T0<=t<Tstar.
%              The last coumn is the value for t>=Tstar.
%        simul = Equal 1 if a stochastic simulation up to period Tmax; 0 if
%              no stoch simulation.
%    simul_mom = If >0 it simulates #simul_mom stochastic paths of size Tmax, 
%              computing for each of them the mean and variance of the path, 
%              and then taking average across simulations. This generates
%              simulated moments from the non-linear solution.
%          irf = Equal 1 if Computation of impulse responses to each shock
%              in the model, using the time varying parameters. It assumes
%              that before the shock the model was at the steady state of
%              the initial parameter value, and that at the same time the
%              stochastic shock hits, the announcement of parameter change 
%              occurs. The resulting IRF is of horizon Tmax-T0+1. Equal 0
%              if no IRF.
% irf_announce = Equal 1 if Computation deterministic IRF with the effect
%              of the "annoncement" of parameter changes at T0, starting
%              from the steady state of the initial parameter value, with
%              horizon Tmax-T0+1; 0 if no IRF.
% irf_announce_simul= If >0 it computes the analogous of irf_announce, but
%              not as a deterministic path but as the avergae across
%              stochastic simulations. Th number of simulations is #irf_announce_simul

% NOTES:
% (1) It assumes that a way to compute the steady state was provided

global M_ oo_ options_
 
M_=M_orig;
oo_=oo_orig;
options_=options_orig;

% Preliminaries
% Position of parameters in M_.param_names
pos_para=strmatch(opt.names_change,M_.param_names,'exact');

%%%%%%%%%%%%%%%
% Solving the model with initial parameters
%%%%%%%%%%%%%%
% Renaming structures
% M_0=M_orig;
% oo_0=oo_orig;
% options_0=options_orig;
% specifying initial parameter
M_.params(pos_para,1)=opt.values_theta(:,1);    
exo_ss = [oo_orig.exo_steady_state; oo_orig.exo_det_steady_state];
    
% Computing steady state
if exist([M_.fname '_steadystate2.m'],'file')==2
    eval(['[ys_0, params0, info] = ' M_.fname '_steadystate2([], exo_ss,M_.params);']);
elseif exist([M_.fname '_steadystate.m'],'file')==2
    eval(['[ys_0, params0, info] = ' M_.fname '_steadystate([], exo_ss,M_,options_);']);
end
% Storing steady state and updated params
oo_.steady_state= ys_0';
M_.params=params0;
% Steady state in the initial parameter values
ss0=oo_.dr.ys;

% If T0>1, also get the solution of the model
% if opt.T0>1
    [dr0,info,M_,oo_] = resol(0,M_,options_,oo_);
    % Transforming the results to the desired format
    [K0,Q0,G0]=sol_dyn_declaration(oo_,M_);
% end


%%%%%%%%%%%%%%%%%
% Solving the model with final parameters
%%%%%%%%%%%%%%%%%
% Renaming structures
M_=M_orig;
oo_=oo_orig;
options_=options_orig;
% specifying initial parameter
M_.params(pos_para,:)=opt.values_theta(:,end);    
% Computing steady state    
if exist([M_.fname '_steadystate2.m'],'file')==2
    eval(['[ys_star, params_star, info] = ' M_.fname '_steadystate2([], exo_ss,M.params);']);
elseif exist([M_.fname '_steadystate.m'],'file')==2
    eval(['[ys_star, params_star, info] = ' M_.fname '_steadystate([], exo_ss, M_,options_);']);
end
% Storing steady state and updated params
oo_.steady_state= ys_star';
M_.params=params_star;

% Get the solution of the model
[dr_star,info,M_,oo_] = resol(0,M_,options_,oo_);
% Transforming the results to the desired format
[Kstar,Qstar,Gstar]=sol_dyn_declaration(oo_,M_);

%%%%%%%%%%%%%%%%%%%%%%%%
% Getting the equilibrium matrices for the transition period
%%%%%%%%%%%%%%%%%%%%%%

% Number of transtion periods
N_trans=opt.Tstar-opt.T0;

if size(opt.values_theta,2)==3 % in the case THETA is always the same in the trasition
    % Renaming structures
    M_=M_orig;
    oo_=oo_orig;
    % Specifying initial parameter
    M_.params(pos_para,1)=opt.values_theta(:,2);  
    % Computing steady state
    if exist([M_.fname '_steadystate2.m'],'file')==2
        eval(['[ys_t, paramst, info] = ' M_.fname '_steadystate2([], exo_ss,M_.params);']);
    elseif exist([M_.fname '_steadystate.m'],'file')==2
        eval(['[ys_t, paramst, info] = ' M_.fname '_steadystate([], exo_ss, M_, options_);']);
    end
    % Storing steady state and updated params
    oo_.steady_state= ys_t';
    M_.params=paramst;
    % Obtaining the matrices with the linearizes equilibrium
    [eq_lin_t]=get_lin_eq(M_,oo_);
	% Storing results in the format used by sol_trans below
    Atrans=repmat(eq_lin_t.A,1,N_trans);
    Btrans=repmat(eq_lin_t.B,1,N_trans);
    Ctrans=repmat(eq_lin_t.C,1,N_trans);
    Dtrans=repmat(eq_lin_t.D,1,N_trans);
    Gammatrans=repmat(eq_lin_t.Gamma,1,N_trans);
else
    if size(opt.values_theta(:,2:end-1),2)~=N_trans;
        error('The number of transtion periods does not coincide with the number of paramter values assigned for that transtion');
    end
    % Creating matrices to put results
    Atrans=zeros(M_.endo_nbr,M_.endo_nbr*N_trans);
    Btrans=zeros(M_.endo_nbr,M_.endo_nbr*N_trans);
    Ctrans=zeros(M_.endo_nbr,M_.endo_nbr*N_trans);
    Dtrans=zeros(M_.endo_nbr,M_.exo_nbr*N_trans);
    Gammatrans=zeros(M_.endo_nbr,N_trans);
    % For each paramter value in the transition
    for t=1:N_trans
        % Renaming structures
        M_=M_orig;
        oo_=oo_orig;
        % Specifying initial parameter
        M_.params(pos_para,:)=opt.values_theta(:,1+t);    
        % Computing steady state
        if exist([M_.fname '_steadystate2.m'],'file')==2
            eval(['[ys_t, paramst, info] = ' M_.fname '_steadystate2([], exo_ss,M_.params);']);
        elseif exist([M_.fname '_steadystate.m'],'file')==2
            eval(['[ys_t, info, paramst] = ' M_.fname '_steadystate([], exo_ss);']);
        end
        % Storing steady state and updated params
        oo_.steady_state= ys_t';
        M_.params=paramst;
        % Obtaining the matrices with the linearizes equilibrium
        [eq_lin_t]=get_lin_eq(M_,oo_);
        % Storing results in the format used by sol_trans below
        Atrans(:,M_.endo_nbr*(t-1)+1:M_.endo_nbr*t)=eq_lin_t.A;
        Btrans(:,M_.endo_nbr*(t-1)+1:M_.endo_nbr*t)=eq_lin_t.B;
        Ctrans(:,M_.endo_nbr*(t-1)+1:M_.endo_nbr*t)=eq_lin_t.C;
        Dtrans(:,M_.exo_nbr*(t-1)+1:M_.exo_nbr*t)=eq_lin_t.D;
        Gammatrans(:,t)=eq_lin_t.Gamma;
    end 
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Finding the solution
%%%%%%%%%%%%%%%%%%%%%%%%%%%
[Ktrans,Qtrans,Gtrans]=sol_trans(N_trans,Atrans,Btrans,Ctrans,Dtrans,Gammatrans,Kstar,Qstar,Gstar);

%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Setting up the results
%%%%%%%%%%%%%%%%%%%%%%%%%%
if opt.T0>1
    Ks=[repmat(K0,1,opt.T0-1) Ktrans repmat(Kstar,1,opt.Tmax-opt.Tstar+1)];
    Qs=[repmat(Q0,1,opt.T0-1) Qtrans repmat(Qstar,1,opt.Tmax-opt.Tstar+1)];
    Gs=[repmat(G0,1,opt.T0-1) Gtrans repmat(Gstar,1,opt.Tmax-opt.Tstar+1)];
else
    Ks=[K0 Ktrans repmat(Kstar,1,opt.Tmax-opt.Tstar+1)];
    Qs=[Q0 Qtrans repmat(Qstar,1,opt.Tmax-opt.Tstar+1)];
    Gs=[G0 Gtrans repmat(Gstar,1,opt.Tmax-opt.Tstar+1)];
end

%%%%%%%%%%%%%%%%
% String results in Dynare structure
%%%%%%%%%%%%%%%%
oo_.sol_para_t.Ks=Ks;
oo_.sol_para_t.Qs=Qs;
oo_.sol_para_t.Gs=Gs;
oo_.sol_para_t.ss0=ss0;
M_.para_t=opt;

%%%%%%%%%%%%%%%%
% Computing simulations
%%%%%%%%%%%%%%%%

% Simulate one stochastic equilibrium path with a random draw from structural shocks
if isfield(opt,'simul')
    if opt.simul ==1 % if stochastic simulation up to period Tmax; 
    exo_simul=mvnrnd(oo_orig.exo_steady_state,M_.Sigma_e,opt.Tmax)';
    y_ini=ss0;
    [y_sim]=simul_model(Ks,Qs,Gs,y_ini,exo_simul,0);
    oo_.sol_para_t.endo_simul=y_sim(1:M_.orig_endo_nbr,:);
    end
end

% Computing simulated moments (mean and variance
if isfield(opt,'simul_mom')
    if opt.simul_mom >0 % if stochastic simulation up to period Tmax; 
        mean_all=nan(M_.orig_endo_nbr,opt.simul_mom);
        std_all =nan(M_.orig_endo_nbr,opt.simul_mom);
        if isfield(opt,'irf_announce_simul')
            if opt.irf_announce_simul >0 % if stochastic simulation up to period Tmax; 
                simul_all_sum=nan(M_.orig_endo_nbr,opt.Tmax+1);
            end
        end
        auxi_ones=kron(ones(1,opt.Tmax+1),ss0);
        for ii=1:max(opt.simul_mom,opt.irf_announce_simul)
            exo_simul_ii=mvnrnd(oo_orig.exo_steady_state,M_.Sigma_e,opt.Tmax)';
            % simulating model with parameter changes
            if isfield(opt,'irf_announce_simul')
                if opt.irf_announce_simul >0 && ii<=opt.irf_announce_simul % if stochastic simulation up to period Tmax; 
                    [y_sim_ii,y_sim_ii_orig]=simul_model(Ks,Qs,Gs,ss0,exo_simul_ii,1);
                end
            else
                [y_sim_ii]=simul_model(Ks,Qs,Gs,ss0,exo_simul_ii,0);
            end
                    
            if ii<=opt.simul_mom
                mean_all(:,ii)=mean(y_sim_ii(:,2:end),2);
                std_all(:,ii)=var(y_sim_ii(:,2:end),0,2).^0.5;
            end
            if isfield(opt,'irf_announce_simul')
                if opt.irf_announce_simul >0 % if stochastic simulation up to period Tmax; 
                    if ii<=opt.irf_announce_simul
                        % Subtracting the steady state if required
                        if opt.irf_announce_type==1
                            y_sim_ii=y_sim_ii-auxi_ones;
                            y_sim_ii_orig=y_sim_ii_orig-auxi_ones;
                        end
                        % Summing responses across smulation to later obtain average
                        if ii==1
                            simul_all_sum=y_sim_ii;
                        else
                            simul_all_sum=simul_all_sum+y_sim_ii;
                        end
                        % Substracting simulation if no parameter change,
                        % in deviation relative to steady state, to keep
                        % the efect of the potential ss change
                        simul_all_sum=simul_all_sum-(y_sim_ii_orig-auxi_ones);
                    end
                end
            end
        end
        oo_.sol_para_t.mean_simul=mean(mean_all,2);
        oo_.sol_para_t.stdev_simul=mean(std_all,2);
        if isfield(opt,'irf_announce_simul')
            if opt.irf_announce_simul >0 % if stochastic simulation up to period Tmax; 
%                 mean_simul=mean(simul_all,3);
%                 for i=1:M_.orig_endo_nbr
%                     eval(['oo_.sol_para_t.irfs.' deblank(M_.endo_names(i,:)) '_announce_simul= mean_simul(i,:) ;']);
%                 end
                mean_simul=simul_all_sum/opt.irf_announce_simul;
                for i=1:M_.orig_endo_nbr
%                     mean_simul_i=mean(squeeze(simul_all_sum(i,:,:)),2);
                    eval(['oo_.sol_para_t.irfs.' deblank(M_.endo_names(i,:)) '_announce_simul= mean_simul(i,:);']);
                end
            end
        end
    end
end

% Impulse responses to structural shocks
if isfield(opt,'irfs')
    if opt.irfs ==1 % if Computation of impulse responses to each shock in the model
    y_ini=ss0;
        for j=1:M_.exo_nbr
            exo_simul_j=zeros(M_.exo_nbr,opt.Tmax);
            exo_simul_j(j,opt.T0)=M_.Sigma_e(j,j)^(0.5);
            [y_sim_j]=simul_model(Ks,Qs,Gs,y_ini,exo_simul_j,0);
            y_sim_j=y_sim_j-kron(ones(1,opt.Tmax+1),y_ini);
            y_sim_j=y_sim_j(:,2:end);
            for i=1:M_.orig_endo_nbr
                eval(['oo_.sol_para_t.irfs.' M_.endo_names{i,:} '_' M_.exo_names{j,:} '= y_sim_j(i,2:end) ;']);
            end
        end
    end
end

% Computing deterministic IRF to the effect of the "annoncement" of parameter changes at T0
if isfield(opt,'irf_announce')
    if opt.irf_announce==1 
        exo_simul=zeros(M_.exo_nbr,opt.Tmax);
        y_ini=ss0;
        [y_sim]=simul_model(Ks,Qs,Gs,y_ini,exo_simul,0);
        if opt.irf_announce_type==1
            y_sim=y_sim-kron(ones(1,opt.Tmax+1),y_ini);
        end
        for i=1:M_.orig_endo_nbr
            eval(['oo_.sol_para_t.irfs.' M_.endo_names{i,:} '_announce= y_sim(i,:) ;']);
        end
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Auxiliary Functions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [K,Q,G]=sol_dyn_declaration(oo_,M_)
% This function takes the model's solution computed by Dynare and sets it
% in the format required by the toolkit (matrices K, Q, G)

% Matrix Q
% Re-sorting the matrix ghx in declaration order
aux_ghx=[oo_.dr.order_var oo_.dr.ghx] ;
aux_ghx=sortrows(aux_ghx);
aux_ghx=aux_ghx(:,2:end);
% creating matrix Q
Q=zeros(M_.endo_nbr,M_.endo_nbr);
Q(:,oo_.dr.state_var)=aux_ghx;
% Matrix G
% Re-sorting the matrix ghu in declaration order
aux_ghu=[oo_.dr.order_var oo_.dr.ghu] ;
aux_ghu=sortrows(aux_ghu);
aux_ghu=aux_ghu(:,2:end);
% creating matrix G
G=aux_ghu;
% Vector K
K=(eye(M_.endo_nbr)-Q)*oo_.dr.ys;

function [Ktrans,Qtrans,Gtrans]=sol_trans(N_trans,Atrans,Btrans,Ctrans,Dtrans,Gammatrans,Kstar,Qstar,Gstar)
% This function computes the solution for the transtion period using 
% backward induction

% Number of variables
[n_y,n_e]=size(Gstar);
% Matrices to store results
Ktrans=nan(n_y,N_trans+1);
Qtrans=nan(n_y,n_y*(N_trans+1));
if n_e>0
    Gtrans=nan(n_y,n_e*(N_trans+1));
else
    Gtrans=[];
end
% Last value of the matrices (erased latter)
Ktrans(:,end)=Kstar;
Qtrans(:,n_y*N_trans+1:end)=Qstar;
if n_e>0
    Gtrans(:,n_e*N_trans+1:end)=Gstar;
end 
% Obtaining the solution by backward induction
for t=N_trans:-1:1
    auxi_t=-sparse(Atrans(:,n_y*(t-1)+1:n_y*t)+Btrans(:,n_y*(t-1)+1:n_y*t)*Qtrans(:,n_y*t+1:n_y*(t+1)))\eye(n_y);
    Qtrans(:,n_y*(t-1)+1:n_y*t)=auxi_t*Ctrans(:,n_y*(t-1)+1:n_y*t);
    Ktrans(:,t)=auxi_t*(Gammatrans(:,t)+Btrans(:,n_y*(t-1)+1:n_y*t)*Ktrans(:,t+1));
    if n_e>0
        Gtrans(:,n_e*(t-1)+1:n_e*t)=auxi_t*Dtrans(:,n_e*(t-1)+1:n_e*t);
    end
end
    
% Cleaning up, to keep only the transition matrices (not the terminal
% point)
Ktrans=Ktrans(:,1:end-1);
Qtrans=Qtrans(:,1:n_y*(N_trans));
if n_e>0
    Gtrans=Gtrans(:,1:n_e*(N_trans));
end


function [eq_lin]=get_lin_eq(M_,oo_)
% This function saves the linearized equilibirum conditions evaluated at 
% the parameter values generated by Dynare, in teh toolkit's format. 
% Let Y be the vector of all endoegnous variebles and U be de vector of
% exogenous variables. The linearized equilibrium condition take the form:
%
%  A*(Y_{t}-Y_ss)+B0*(Y_{t+1}-Y_ss)+C*(Y_{t-1}-Y_ss)+D*U_{t}
%
% or,
%
%  A*Y_{t}+B*Y_{t+1}+C*Y_{t-1}+D*U_{t}+Gamma, with Gamma=-(A+B+C)*Y_ss,
%
% where Y_ss is the steady state of Y_{t}, and the matrices A0, A1, B0, D0 
% are the Jacobian matrices with the derivatives of F with respect to,
% respectively, Y_{t+1}, Y_{t}, Y_{t-1}, and U_{t}, evaluated at the state
% state. It is assumed that U_{t} equals 0 in steady state.
%
% The inputs are the structures M_ and oo_ produced by Dynare. The output
% is a structure with with content:
%     eq_lin.A, eq_lin.B, eq_lin.C, eq_lin.D, eq_lin.Gamma
% representing, respectively, the matrices A, B, C, D and Gamma.
%
% NOTES: 
% (1) This function assume that the model has at most one lead and one
% lag of the variables in Y
% (2) The function was created following the m-file 'stochastic_solvers.m'
% the in included in the 'matlab' folder od Dynare.

klen1 = M_.maximum_lag + M_.maximum_lead + 1;
exo_simul1 = [repmat(oo_.exo_steady_state',klen1,1) repmat(oo_.exo_det_steady_state',klen1,1)];
iyv1 = M_.lead_lag_incidence';
iyv1 = iyv1(:);
iyr01 = find(iyv1) ;
it_1 = M_.maximum_lag + 1 ;

if M_.exo_nbr == 0
    oo_.exo_steady_state = [] ;
end

it_1 = M_.maximum_lag + 1;
% keyboard;
z1 = repmat(oo_.steady_state,1,klen1);

[~,jaco_] = feval([M_.fname '.dynamic'],z1(iyr01),exo_simul1,M_.params, oo_.steady_state, it_1);

% Matrix with Y_{t}
eq_lin.A=zeros(M_.endo_nbr,M_.endo_nbr);
auxi_cont=find(M_.lead_lag_incidence(2,:));
eq_lin.A(:,auxi_cont)=jaco_(:,M_.lead_lag_incidence(2,auxi_cont));

% Matrix with Y_{t+1}
eq_lin.B=zeros(M_.endo_nbr,M_.endo_nbr);
auxi_cont=find(M_.lead_lag_incidence(3,:));
eq_lin.B(:,auxi_cont)=jaco_(:,M_.lead_lag_incidence(3,auxi_cont));

% Matrix with Y_{t-1}
eq_lin.C=zeros(M_.endo_nbr,M_.endo_nbr);
auxi_cont=find(M_.lead_lag_incidence(1,:));
eq_lin.C(:,auxi_cont)=jaco_(:,M_.lead_lag_incidence(1,auxi_cont));

% Matrix with U_{t}
eq_lin.D=jaco_(:,end-M_.exo_nbr+1:end);

% Matrix with Constant
% eq_lin.Gamma=-(eq_lin.A+eq_lin.B+eq_lin.C)*oo_.dr.ys;
eq_lin.Gamma=-(eq_lin.A+eq_lin.B+eq_lin.C)*oo_.steady_state';

function [y_sim,y_sim_ini]=simul_model(Ks,Qs,Gs,y_ini,e_sim,do_simul_ini,elim_init)
% This function simulates the model using the sequence of solutions in 
% Ks,Qs,Gs, starting from a point y_ini, asuming a sequence for teh
% exogenous variables given by e_sim
[n_y,T]=size(Ks);
n_e=size(Gs,2)/T;

y_sim=nan(n_y,T+1);

y_sim(:,1)=y_ini;

if exist('do_simul_ini')==1 && nargout>1
    if do_simul_ini==1
        exist_ini=1;
        y_sim_ini=y_sim;
    end
else
    exist_ini=0;
end

for t=1:size(e_sim,2) % BE CAREFUl: For the matrices Ks, Qs, Gs, t is actually teh time period, but for y_sim is t+1 (due to the initial point)
    if T>1
        y_sim(:,t+1)=Ks(:,t)+Qs(:,(t-1)*n_y+1:t*n_y)*y_sim(:,t);    
        if n_e>0
            y_sim(:,t+1)=y_sim(:,t+1)+Gs(:,(t-1)*n_e+1:t*n_e)*e_sim(:,t);
        end
        if exist_ini==1
            y_sim_ini(:,t+1)=Ks(:,1)+Qs(:,1:n_y)*y_sim_ini(:,t);    
            if n_e>0
                y_sim_ini(:,t+1)=y_sim_ini(:,t+1)+Gs(:,1:n_e)*e_sim(:,t);
            end
        end
    else
        y_sim(:,t+1)=Ks+Qs*y_sim(:,t);    
        if n_e>0
            y_sim(:,t+1)=y_sim(:,t+1)+Gs*e_sim(:,t);
        end
        if exist_ini==1
            y_sim_ini(:,t+1)=y_sim(:,t+1);
        end
    end
end

%eliminting initial point
if nargin>6
    if elim_init==1
        y_sim=y_sim(:,2:end);
        if exist_ini==1
            y_sim_ini=y_sim;
        end
    end
end