dr1_PI.m 19.2 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
function [dr,info,M_,options_,oo_] = dr1_PI(dr,task,M_,options_,oo_)
% function [dr,info,M_,options_,oo_] = dr1_PI(dr,task,M_,options_,oo_)
% Computes the reduced form solution of a rational expectation model first
% order
% approximation of the Partial Information stochastic model solver around the deterministic steady state).
% Prepares System as
%        A0*E_t[y(t+1])+A1*y(t)=A2*y(t-1)+c+psi*eps(t)
% with z an exogenous variable process.
% and calls PI_Gensys.m solver 
% based on Pearlman et al 1986 paper and derived from
% C.Sims' gensys linear solver.
% to return solution in format
%       [s(t)' x(t)' E_t x(t+1)']'=G1pi [s(t-1)' x(t-1)' x(t)]'+C+impact*eps(t),
%
% INPUTS
%   dr         [matlab structure] Decision rules for stochastic simulations.
%   task       [integer]          if task = 0 then dr1 computes decision rules.
%                                 if task = 1 then dr1 computes eigenvalues.
%   M_         [matlab structure] Definition of the model.           
%   options_   [matlab structure] Global options.
%   oo_        [matlab structure] Results 
%    
% OUTPUTS
%   dr         [matlab structure] Decision rules for stochastic simulations.
%   info       [integer]          info=1: the model doesn't define current variables uniquely
%                                 info=2: problem in mjdgges.dll info(2) contains error code. 
%                                 info=3: BK order condition not satisfied info(2) contains "distance"
%                                         absence of stable trajectory.
%                                 info=4: BK order condition not satisfied info(2) contains "distance"
%                                         indeterminacy.
%                                 info=5: BK rank condition not satisfied.
%   M_         [matlab structure]            
%   options_   [matlab structure]
%   oo_        [matlab structure]
%  
% ALGORITHM
%   ...
%    
% SPECIAL REQUIREMENTS
%   none.
%  

43
% Copyright (C) 1996-2012 Dynare Team
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
%
% 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 <http://www.gnu.org/licenses/>.

60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
global lq_instruments;
info = 0;

options_ = set_default_option(options_,'qz_criterium',1.000001);

xlen = M_.maximum_endo_lead + M_.maximum_endo_lag + 1;

if (options_.aim_solver == 1)
    options_.aim_solver == 0;
    warning('You can not use AIM with Part Info solver. AIM ignored'); 
end
if (options_.order > 1) 
    warning('You can not use order higher than 1 with Part Info solver. Order 1 assumed');
    options_.order =1;
end
75

76
% expanding system for Optimal Linear Regulator
77
if options_.ramsey_policy && options_.ACES_solver == 0
78
79
80
81
82
83
84
85
86
    if isfield(M_,'orig_model')
        orig_model = M_.orig_model;
        M_.endo_nbr = orig_model.endo_nbr;
        M_.endo_names = orig_model.endo_names;
        M_.lead_lag_incidence = orig_model.lead_lag_incidence;
        M_.maximum_lead = orig_model.maximum_lead;
        M_.maximum_endo_lead = orig_model.maximum_endo_lead;
        M_.maximum_lag = orig_model.maximum_lag;
        M_.maximum_endo_lag = orig_model.maximum_endo_lag;
87
    end
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
    old_solve_algo = options_.solve_algo;
    %  options_.solve_algo = 1;
    oo_.steady_state = dynare_solve('ramsey_static',oo_.steady_state,0,M_,options_,oo_,it_);
    options_.solve_algo = old_solve_algo;
    [junk,junk,multbar] = ramsey_static(oo_.steady_state,M_,options_,oo_,it_);
    [jacobia_,M_] = ramsey_dynamic(oo_.steady_state,multbar,M_,options_,oo_,it_);
    klen = M_.maximum_lag + M_.maximum_lead + 1;
    dr.ys = [oo_.steady_state;zeros(M_.exo_nbr,1);multbar];
else
    klen = M_.maximum_lag + M_.maximum_lead + 1;
    iyv = M_.lead_lag_incidence';
    iyv = iyv(:);
    iyr0 = find(iyv) ;
    it_ = M_.maximum_lag + 1 ;
    
    if M_.exo_nbr == 0
        oo_.exo_steady_state = [] ;
105
    end
106
    
107

108
109
110
111
    if options_.ACES_solver == 1
        sim_ruleids=[];
        tct_ruleids=[];
        if  size(M_.equations_tags,1)>0  % there are tagged equations, check if they are aceslq rules
112
            for teq=1:size(M_.equations_tags,1)
113
114
115
116
117
118
                if strcmp(M_.equations_tags(teq,3),'aceslq_sim_rule')
                    sim_ruleids=[sim_ruleids cell2mat(M_.equations_tags(teq,1))]
                end
                if strcmp(M_.equations_tags(teq,3),'aceslq_tct_rule')
                    tct_ruleids=[tct_ruleids cell2mat(M_.equations_tags(teq,1))]
                end
119
120
            end
        end
121
122
123
124
125
126
127
128
129
130
131
        lq_instruments.sim_ruleids=sim_ruleids;
        lq_instruments.tct_ruleids=tct_ruleids;
        %if isfield(lq_instruments,'xsopt_SS') %% changed by BY
        [junk, lq_instruments.xsopt_SS,lq_instruments.lmopt_SS,s2,check] = opt_steady_get;%% changed by BY
        [qc, DYN_Q] = QPsolve(lq_instruments, s2, check); %% added by BY
        z = repmat(lq_instruments.xsopt_SS,1,klen);
    else
        z = repmat(dr.ys,1,klen);
    end
    z = z(iyr0) ;
    [junk,jacobia_] = feval([M_.fname '_dynamic'],z,[oo_.exo_simul ...
132
                        oo_.exo_det_simul], M_.params, dr.ys, it_);
133

134
135
    if options_.ACES_solver==1 && (length(sim_ruleids)>0 || length(tct_ruleids)>0 )
        if length(sim_ruleids)>0
136
137
138
            sim_rule=jacobia_(sim_ruleids,:);
            % uses the subdirectory - BY
            save ([ACES_DirectoryName,'/',M_.fname '_ACESLQ_sim_rule.txt'], 'sim_rule', '-ascii', '-double', '-tabs');
139
140
        end
        if length(tct_ruleids)>0
141
142
143
144
            tct_rule=jacobia_(tct_ruleids,:);
            % uses the subdirectory - BY
            save ([ACES_DirectoryName,'/',M_.fname '_ACESLQ_tct_rule.txt'], 'tct_rule', '-ascii', '-double', '-tabs');
        end
145
146
147
148
149
        aces_ruleids=union(tct_ruleids,sim_ruleids);
        j_size=size(jacobia_,1);
        j_rows=1:j_size;
        j_rows = setxor(j_rows,aces_ruleids);
        jacobia_=jacobia_(j_rows ,:);
150
    end
151

152
153
154
155
156
157
end

if options_.debug
    save([M_.fname '_debug.mat'],'jacobia_')
end

158
dr=set_state_space(dr,M_,options_);
159
160
161
162
163
164
165
166
kstate = dr.kstate;
nstatic = dr.nstatic;
nfwrd = dr.nfwrd;
npred = dr.npred;
nboth = dr.nboth;
order_var = dr.order_var;
nd = size(kstate,1);
nz = nnz(M_.lead_lag_incidence);
167

168
169
170
171
172
173
174
175
sdyn = M_.endo_nbr - nstatic;

k0 = M_.lead_lag_incidence(M_.maximum_endo_lag+1,order_var);
k1 = M_.lead_lag_incidence(find([1:klen] ~= M_.maximum_endo_lag+1),:);

if (options_.aim_solver == 1) 
    error('Anderson and Moore AIM solver is not compatible with Partial Information models');
end % end if useAIM and...
176
177
178
179
180

    % If required, try PCL86 solver, that is, if not the check being
    % performed only and if it is 1st order 
    % create sparse, extended jacobia AA:
    nendo=M_.endo_nbr; % = size(aa,1)
181
    
182

183
    if(options_.ACES_solver==1)
184
        %if ~isfield(lq_instruments,'names')
185
        if isfield(options_,'instruments')
186
            lq_instruments.names=options_.instruments;
187
        end
188
189
190
        %end
        if isfield(lq_instruments,'names')
            num_inst=size(lq_instruments.names,1);
191
            if ~isfield(lq_instruments,'inst_var_indices') && num_inst>0
192
193
194
195
196
197
198
199
200
201
202
203
204
205
                for i=1:num_inst
                    i_tmp = strmatch(deblank(lq_instruments.names(i,:)),M_.endo_names,'exact');
                    if isempty(i_tmp)
                        error (['One of the specified instrument variables does not exist']) ;
                    else
                        i_var(i) = i_tmp;
                    end
                end
                lq_instruments.inst_var_indices=i_var;
            elseif size(lq_instruments.inst_var_indices)>0
                i_var=lq_instruments.inst_var_indices;
                if ~num_inst
                    num_inst=size(lq_instruments.inst_var_indices);
                end
206
            else
207
208
                i_var=[];
                num_inst=0;
209
            end
210
            if size(i_var,2)>0 && size(i_var,2)==num_inst
211
212
213
214
215
216
217
218
219
220
                m_var=zeros(nendo,1);
                for i=1:nendo
                    if isempty(find(i_var==i))
                        m_var(i)=i;
                    end
                end
                m_var=nonzeros(m_var);
                lq_instruments.m_var=m_var;
            else
                error('WARNING: There are no instrumnets for ACES!');
221
            end
222
223
        else %if(options_.ACES_solver==1)
            error('WARNING: There are no instrumnets for ACES!');
224
225
        end
    end
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240

    % find size xlen of the state vector Y and of A0, A1 and A2 transition matrices:
    % it is the sum the all i variables's lag/lead representations,
    % for each variable i representation being defined as:
    % Max (i_lags-1,0)+ Max (i_leads-1,0)+1
    % so that if variable x appears with 2 lags and 1 lead, and z
    % with 2 lags and 3 leads, the size of the state space is:
    % 1+0+1   +   1+2+1   =6
    % e.g. E_t Y(t+1)=
    %     E_t x(t)
    %     E_t x(t+1)
    %     E_t z(t)
    %     E_t z(t+1)
    %     E_t z(t+2)
    %     E_t z(t+3) 
241
242
    
    % partition jacobian:
243
244
    jlen=dr.nspred+dr.nsfwrd+M_.endo_nbr+M_.exo_nbr; % length of jacobian
    PSI=-jacobia_(:, jlen-M_.exo_nbr+1:end); % exog
245
                                             % first transpose M_.lead_lag_incidence';
246
247
248
    lead_lag=M_.lead_lag_incidence';
    max_lead_lag=zeros(nendo,2); % lead/lag representation in Y for each endogenous variable i
    if ( M_.maximum_lag <= 1) && (M_.maximum_lead <= 1)
249
250
251
252
253
254
255
256
257
258
        xlen=size(jacobia_,1);%nendo;
        AA0=zeros(xlen,xlen);  % empty A0
        AA2=AA0; % empty A2 and A3
        AA3=AA0;
        if xlen==nendo % && M_.maximum_lag <=1 && M_.maximum_lead <=1 % apply a shortcut
            AA1=jacobia_(:,npred+1:npred+nendo); 
            if M_.maximum_lead ==1
                fnd = find(lead_lag(:,M_.maximum_lag+2));
                AA0(:, fnd)= jacobia_(:,nonzeros(lead_lag(:,M_.maximum_lag+2))); %forwd jacobian
            end
259
            if npred>0 && M_.maximum_lag ==1  
260
261
262
263
                fnd = find(lead_lag(:,1));
                AA2(:, fnd)= jacobia_(:,nonzeros(lead_lag(:,1))); %backward
            end
        elseif options_.ACES_solver==1 % more endo vars than equations in jacobia_
264
        if nendo-xlen==num_inst
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
            PSI=[PSI;zeros(num_inst, M_.exo_nbr)];
            % AA1 contemporary
            AA_all=jacobia_(:,npred+1:npred+nendo); 
            AA1=AA_all(:,lq_instruments.m_var); % endo without instruments
            lq_instruments.ij1=AA_all(:,lq_instruments.inst_var_indices); %  instruments only
            lq_instruments.B1=-[lq_instruments.ij1; eye(num_inst)];
            AA1=[AA1, zeros(xlen,num_inst); zeros(num_inst,xlen), eye(num_inst)];
            %PSI=[PSI; zeros(num_inst,M_.exo_nbr)];
            if M_.maximum_lead ==1 % AA0 forward looking
                AA_all(:,:)=0.0;
                fnd = find(lead_lag(:,M_.maximum_lag+2));
                AA_all(:, fnd)= jacobia_(:,nonzeros(lead_lag(:,M_.maximum_lag+2))); %forwd jacobian
                AA0=AA_all(:,lq_instruments.m_var);
                lq_instruments.ij0=AA_all(:,lq_instruments.inst_var_indices); %  instruments only
                lq_instruments.B0=[lq_instruments.ij0; eye(num_inst)];
                AA0=[AA0, zeros(xlen,num_inst); zeros(num_inst,xlen+num_inst)];
            end
282
            if npred>0 && M_.maximum_lag ==1  
283
284
285
286
287
288
289
290
291
292
293
                AA_all(:,:)=0.0;
                fnd = find(lead_lag(:,1));
                AA_all(:, fnd)= jacobia_(:,nonzeros(lead_lag(:,1))); %backward
                AA2=AA_all(:,lq_instruments.m_var);
                lq_instruments.ij2=AA_all(:,lq_instruments.inst_var_indices); %  instruments only
                lq_instruments.B2=[lq_instruments.ij2; eye(num_inst)];
                AA2=[AA2, lq_instruments.ij2 ; zeros(num_inst,xlen+num_inst)];
            end
        else
            error('ACES number of instruments does match');
        end
294
        else
295
            error('More than one lead or lag in the jabian');
296
        end
297
298
299
300
301
302
303
304
305
        if M_.orig_endo_nbr<nendo
            % findif there are any expecatations at time t
            exp_0= strmatch('AUX_EXPECT_LEAD_0_', M_.endo_names);
            num_exp_0=size(exp_0,1);
            if num_exp_0>0
                AA3(:,exp_0)=AA1(:,exp_0);
                XX0=zeros(nendo,num_exp_0);
                AA1(:,exp_0)=XX0(:,[1:num_exp_0])
            end
306
307
        end
    end
308
309
310
311
312
    PSI=-[[zeros(xlen-nendo,M_.exo_nbr)];[jacobia_(:, jlen-M_.exo_nbr+1:end)]]; % exog
    cc=0;
    NX=M_.exo_nbr; % no of exogenous varexo shock variables.
    NETA=nfwrd+nboth; % total no of exp. errors  set to no of forward looking equations
    FL_RANK=rank(AA0); % nfwrd+nboth; % min total no of forward looking equations and vars
313
    
314
    try
315
316
317
318
319
320
321
322
323
324
        % call [G1pi,C,impact,nmat,TT1,TT2,gev,eu]=PI_gensys(a0,a1,a2,c,PSI,NX,NETA,NO_FL_EQS)
        % System given as
        %        a0*E_t[y(t+1])+a1*y(t)=a2*y(t-1)+c+psi*eps(t)
        % with eps an exogenous variable process.
        % Returned system is
        %       [s(t)' x(t)' E_t x(t+1)']'=G1pi [s(t-1)' x(t-1)' x(t)]'+C+impact*eps(t),
        %  and (a) the matrix nmat satisfying   nmat*E_t z(t)+ E_t x(t+1)=0
        %      (b) matrices TT1, TT2  that relate y(t) to these states:
        %      y(t)=[TT1 TT2][s(t)' x(t)']'.

325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
        if(options_.ACES_solver==1)
            if isfield(lq_instruments,'xsopt_SS')
                SSbar= diag([lq_instruments.xsopt_SS(m_var)]);% lq_instruments.xsopt_SS(lq_instruments.inst_var_indices)]);
                insSSbar=repmat(lq_instruments.xsopt_SS(lq_instruments.inst_var_indices)',nendo-num_inst,1);
            else
                SSbar= diag([dr.ys(m_var)]);%; dr.ys(lq_instruments.inst_var_indices)]);%(oo_.steady_state);
                insSSbar=repmat(dr.ys(lq_instruments.inst_var_indices)',nendo-num_inst,1);
            end
            SSbar=diag([diag(SSbar);diag(eye(num_inst))]);
            insSSbar=[insSSbar;diag(eye(num_inst))];

            AA0=AA0*SSbar;
            AA1=AA1*SSbar;
            AA2=AA2*SSbar;
            lq_instruments.B1=(lq_instruments.B1).*insSSbar;
        end
        %% for expectational models when complete
        if any(AA3)
            AA3=AA3*SSbar;
            [G1pi,CC,impact,nmat,TT1,TT2,gev,eu, DD, E2,E5, GAMMA, FL_RANK]=PI_gensysEXP(AA0,AA1,-AA2,AA3,cc,PSI,NX,NETA,FL_RANK, M_, options_);
345
        else
346
            [G1pi,CC,impact,nmat,TT1,TT2,gev,eu, DD, E2,E5, GAMMA, FL_RANK]=PI_gensys(AA0,AA1,-AA2,AA3,cc,PSI,NX,NETA,FL_RANK, M_, options_);
347
        end
348

349
        % reuse some of the bypassed code and tests that may be needed 
Sébastien Villemot's avatar
Sébastien Villemot committed
350
        if (eu(1) ~= 1 || eu(2) ~= 1) && options_.ACES_solver==0
351
352
353
354
            info(1) = abs(eu(1)+eu(2));
            info(2) = 1.0e+8;
            %     return
        end
355
        
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
        dr.PI_ghx=G1pi;
        dr.PI_ghu=impact;
        dr.PI_TT1=TT1;
        dr.PI_TT2=TT2;
        dr.PI_nmat=nmat;
        dr.PI_CC=CC;
        dr.PI_gev=gev;
        dr.PI_eu=eu;
        dr.PI_FL_RANK=FL_RANK;
        %dr.ys=zeros(nendo); % zero steady state
        dr.ghx=G1pi;
        dr.ghu=impact;
        dr.eigval = eig(G1pi);
        dr.rank=FL_RANK;
        
        if options_.ACES_solver==1
            betap=options_.planner_discount;
            sigma_cov=M_.Sigma_e;
            % get W - BY
            W=(1-betap)*GAMMA'*DYN_Q*GAMMA;
            %W=[0]
            ACES.A=G1pi;
            ACES.C=impact; % (:,1);
            ACES.D=DD; %=impact (:,20);
            ACES.E2=E2;
            ACES.E5=E5;
            ACES.GAMMA=GAMMA;
            ACES_M=size(G1pi,2)-FL_RANK;
            ACES_NM=FL_RANK;
            ACES.M=ACES_M;
            ACES.NM=FL_RANK;
            % added by BY
            ACES.Q=DYN_Q;
            ACES.W=W;
            NY=nendo-num_inst;
            
            % save the followings in a subdirectory - BY    
            save ([ACES_DirectoryName,'/',M_.fname '_ACESLQ_Matrices'], 'ACES');
            save ([ACES_DirectoryName,'/',M_.fname '_ACESLQ_GAMMA'], 'GAMMA');
            save ([ACES_DirectoryName,'/',M_.fname '_ACESLQ_A.txt'], 'G1pi', '-ascii', '-double', '-tabs');
            save ([ACES_DirectoryName,'/',M_.fname '_ACESLQ_C.txt'], 'impact','-ascii', '-double', '-tabs');
            save ([ACES_DirectoryName,'/',M_.fname '_ACESLQ_D.txt'], 'DD', '-ascii', '-double', '-tabs');
            save ([ACES_DirectoryName,'/',M_.fname '_ACESLQ_E2.txt'], 'E2', '-ascii', '-double', '-tabs');
            save ([ACES_DirectoryName,'/',M_.fname '_ACESLQ_E5.txt'], 'E5', '-ascii', '-double', '-tabs');
            save ([ACES_DirectoryName,'/',M_.fname '_ACESLQ_GAMMA.txt'], 'GAMMA', '-ascii', '-double', '-tabs');
            %save ([M_.fname '_ACESLQ_M.txt'], 'ACES_M', '-ascii', '-tabs');
            %save ([M_.fname '_ACESLQ_NM.txt'], 'ACES_NM', '-ascii', '-tabs');
            %save ([M_.fname '_ACESLQ_betap.txt'], 'betap', '-ascii', '-tabs');
            %save ([M_.fname '_ACESLQ_NI.txt'], 'num_inst', '-ascii', '-tabs');
            %save ([M_.fname '_ACESLQ_ND.txt'], 'NX', '-ascii', '-tabs');
            %save ([M_.fname '_ACESLQ_NY.txt'], 'NY', '-ascii', '-tabs');
            ACES_VARS=M_.endo_names;
            save ([ACES_DirectoryName,'/',M_.fname '_ACESLQ_VARS.txt'], 'ACES_VARS', '-ascii', '-tabs');
            % added by BY
            % save the char array ACES_VARS into .txt as it is
            fid = fopen(strcat(ACES_DirectoryName,'/',M_.fname,'_ACESLQ_VARnames.txt'),'wt');
            ACES_VARS =[ACES_VARS repmat(sprintf('\n'),size(ACES_VARS,1),1)];
            fwrite(fid,ACES_VARS.');
            fclose(fid);
            % save as integers
            fid = fopen(strcat(ACES_DirectoryName,'/',M_.fname,'_ACESLQ_M.txt'),'wt');
            fprintf(fid,'%d\n',ACES_M);
            fclose(fid);
            fid = fopen(strcat(ACES_DirectoryName,'/',M_.fname,'_ACESLQ_NM.txt'),'wt');
            fprintf(fid,'%d\n',ACES_NM);
            fclose(fid);
            fid = fopen(strcat(ACES_DirectoryName,'/',M_.fname,'_ACESLQ_betap.txt'),'wt');
            fprintf(fid,'%d\n',betap);
            fclose(fid);
            fid = fopen(strcat(ACES_DirectoryName,'/',M_.fname,'_ACESLQ_NI.txt'),'wt');
            fprintf(fid,'%d\n',num_inst);
            fclose(fid);
            fid = fopen(strcat(ACES_DirectoryName,'/',M_.fname,'_ACESLQ_ND.txt'),'wt');
            fprintf(fid,'%d\n',NX);
            fclose(fid);
            fid = fopen(strcat(ACES_DirectoryName,'/',M_.fname,'_ACESLQ_NY.txt'),'wt');
            fprintf(fid,'%d\n',NY);
            fclose(fid);
            save ([ACES_DirectoryName,'/',M_.fname '_ACESLQ_Q.txt'], 'DYN_Q', '-ascii', '-double', '-tabs');
            save ([ACES_DirectoryName,'/',M_.fname '_ACESLQ_W.txt'], 'W', '-ascii', '-double', '-tabs');
            save ([ACES_DirectoryName,'/',M_.fname '_ACESLQ_SIGMAE.txt'], 'sigma_cov', '-ascii', '-double', '-tabs');
        end
438
439

    catch
440
441
442
443
444
        lerror=lasterror;
        if options_.ACES_solver==1
            disp('Problem with using Part Info ACES solver:');
            error(lerror.message);
        else
Sébastien Villemot's avatar
Sébastien Villemot committed
445
446
            disp('Problem with using Part Info solver');
            error(lerror.message);
447
        end
448
    end
449
450
451
452

    % TODO: 
    % if options_.loglinear == 1
    % if exogenous deterministic variables
453
454
    
    return;