diff --git a/matlab/get_identification_jacobians.m b/matlab/get_identification_jacobians.m
index 91afeecc5f2af70ffa4756434ad7489da07a7923..0620b7f4805e2963317ebe910e21a093f3767e0f 100644
--- a/matlab/get_identification_jacobians.m
+++ b/matlab/get_identification_jacobians.m
@@ -442,14 +442,16 @@ if ~no_identification_minimal
dMINIMAL = [];
else
% Derive and check minimal state vector of first-order
- [CheckCO,minnx,minA,minB,minC,minD,dminA,dminB,dminC,dminD] = get_minimal_state_representation(oo.dr.ghx(oo.dr.pruned.indx,:),... %A
- oo.dr.ghu(oo.dr.pruned.indx,:),... %B
- oo.dr.ghx(oo.dr.pruned.indy,:),... %C
- oo.dr.ghu(oo.dr.pruned.indy,:),... %D
- oo.dr.derivs.dghx(oo.dr.pruned.indx,:,:),... %dA
- oo.dr.derivs.dghu(oo.dr.pruned.indx,:,:),... %dB
- oo.dr.derivs.dghx(oo.dr.pruned.indy,:,:),... %dC
- oo.dr.derivs.dghu(oo.dr.pruned.indy,:,:)); %dD
+ SYS.A = oo.dr.ghx(oo.dr.pruned.indx,:);
+ SYS.dA = oo.dr.derivs.dghx(oo.dr.pruned.indx,:,:);
+ SYS.B = oo.dr.ghu(oo.dr.pruned.indx,:);
+ SYS.dB = oo.dr.derivs.dghu(oo.dr.pruned.indx,:,:);
+ SYS.C = oo.dr.ghx(oo.dr.pruned.indy,:);
+ SYS.dC = oo.dr.derivs.dghx(oo.dr.pruned.indy,:,:);
+ SYS.D = oo.dr.ghu(oo.dr.pruned.indy,:);
+ SYS.dD = oo.dr.derivs.dghu(oo.dr.pruned.indy,:,:);
+ [CheckCO,minnx,SYS] = get_minimal_state_representation(SYS,1);
+
if CheckCO == 0
warning_KomunjerNg = 'WARNING: Komunjer and Ng (2011) failed:\n';
warning_KomunjerNg = [warning_KomunjerNg ' Conditions for minimality are not fullfilled:\n'];
@@ -457,6 +459,10 @@ if ~no_identification_minimal
fprintf(warning_KomunjerNg); %use sprintf to have line breaks
dMINIMAL = [];
else
+ minA = SYS.A; dminA = SYS.dA;
+ minB = SYS.B; dminB = SYS.dB;
+ minC = SYS.C; dminC = SYS.dC;
+ minD = SYS.D; dminD = SYS.dD;
%reshape into Magnus-Neudecker Jacobians, i.e. dvec(X)/dp
dminA = reshape(dminA,size(dminA,1)*size(dminA,2),size(dminA,3));
dminB = reshape(dminB,size(dminB,1)*size(dminB,2),size(dminB,3));
diff --git a/matlab/get_minimal_state_representation.m b/matlab/get_minimal_state_representation.m
index 4d349eeb5bfbd814c37ce5bbe92997880203a78f..fa25648842ee6a08a2061593bda14a485eb267a8 100644
--- a/matlab/get_minimal_state_representation.m
+++ b/matlab/get_minimal_state_representation.m
@@ -1,40 +1,55 @@
-function [CheckCO,minnx,minA,minB,minC,minD,dminA,dminB,dminC,dminD] = get_minimal_state_representation(A,B,C,D,dA,dB,dC,dD)
-% [CheckCO,minnx,minA,minB,minC,minD,dminA,dminB,dminC,dminD] = get_minimal_state_representation(A,B,C,D,dA,dB,dC,dD)
-% Derives and checks the minimal state representation of the ABCD
-% representation of a state space model
+function [CheckCO,minns,minSYS] = get_minimal_state_representation(SYS, derivs_flag)
+% Derives and checks the minimal state representation
+% Let x = A*x(-1) + B*u and y = C*x(-1) + D*u be a linear state space
+% system, then this function computes the following representation
+% xmin = minA*xmin(-1) + minB*u and and y=minC*xmin(-1) + minD*u
+%
% -------------------------------------------------------------------------
% INPUTS
-% A: [endo_nbr by endo_nbr] Transition matrix from Kalman filter
-% for all endogenous declared variables, in DR order
-% B: [endo_nbr by exo_nbr] Transition matrix from Kalman filter
-% mapping shocks today to endogenous variables today, in DR order
-% C: [obs_nbr by endo_nbr] Measurement matrix from Kalman filter
-% linking control/observable variables to states, in DR order
-% D: [obs_nbr by exo_nbr] Measurement matrix from Kalman filter
-% mapping shocks today to controls/observables today, in DR order
-% dA: [endo_nbr by endo_nbr by totparam_nbr] in DR order
+% SYS [structure]
+% with the following necessary fields:
+% A: [nspred by nspred] in DR order
+% Transition matrix for all state variables
+% B: [nspred by exo_nbr] in DR order
+% Transition matrix mapping shocks today to states today
+% C: [varobs_nbr by nspred] in DR order
+% Measurement matrix linking control/observable variables to states
+% D: [varobs_nbr by exo_nbr] in DR order
+% Measurement matrix mapping shocks today to controls/observables today
+% and optional fields:
+% dA: [nspred by nspred by totparam_nbr] in DR order
% Jacobian (wrt to all parameters) of transition matrix A
-% dB: [endo_nbr by exo_nbr by totparam_nbr] in DR order
+% dB: [nspred by exo_nbr by totparam_nbr] in DR order
% Jacobian (wrt to all parameters) of transition matrix B
-% dC: [obs_nbr by endo_nbr by totparam_nbr] in DR order
+% dC: [varobs_nbr by nspred by totparam_nbr] in DR order
% Jacobian (wrt to all parameters) of measurement matrix C
-% dD: [obs_nbr by exo_nbr by totparam_nbr] in DR order
+% dD: [varobs_nbr by exo_nbr by totparam_nbr] in DR order
% Jacobian (wrt to all parameters) of measurement matrix D
+% derivs_flag [scalar]
+% (optional) indicator whether to output parameter derivatives
% -------------------------------------------------------------------------
% OUTPUTS
-% CheckCO: [scalar] indicator, equals to 1 if minimal state representation is found
-% minnx: [scalar] length of minimal state vector
-% minA: [minnx by minnx] Transition matrix A for evolution of minimal state vector
-% minB: [minnx by exo_nbr] Transition matrix B for evolution of minimal state vector
-% minC: [obs_nbr by minnx] Measurement matrix C for evolution of controls, depending on minimal state vector only
-% minD: [obs_nbr by minnx] Measurement matrix D for evolution of controls, depending on minimal state vector only
-% dminA: [minnx by minnx by totparam_nbr] in DR order
+% CheckCO: [scalar]
+% equals to 1 if minimal state representation is found
+% minns: [scalar]
+% length of minimal state vector
+% SYS [structure]
+% with the following fields:
+% minA: [minns by minns] in DR-order
+% transition matrix A for evolution of minimal state vector
+% minB: [minns by exo_nbr] in DR-order
+% transition matrix B for evolution of minimal state vector
+% minC: [varobs_nbr by minns] in DR-order
+% measurement matrix C for evolution of controls, depending on minimal state vector only
+% minD: [varobs_nbr by minns] in DR-order
+% measurement matrix D for evolution of controls, depending on minimal state vector only
+% dminA: [minns by minns by totparam_nbr] in DR order
% Jacobian (wrt to all parameters) of transition matrix minA
-% dminB: [minnx by exo_nbr by totparam_nbr] in DR order
+% dminB: [minns by exo_nbr by totparam_nbr] in DR order
% Jacobian (wrt to all parameters) of transition matrix minB
-% dminC: [obs_nbr by minnx by totparam_nbr] in DR order
+% dminC: [varobs_nbr by minns by totparam_nbr] in DR order
% Jacobian (wrt to all parameters) of measurement matrix minC
-% dminD: [obs_nbr by exo_nbr by totparam_nbr] in DR order
+% dminD: [varobs_nbr by u_nbr by totparam_nbr] in DR order
% Jacobian (wrt to all parameters) of measurement matrix minD
% -------------------------------------------------------------------------
% This function is called by
@@ -42,8 +57,9 @@ function [CheckCO,minnx,minA,minB,minC,minD,dminA,dminB,dminC,dminD] = get_minim
% -------------------------------------------------------------------------
% This function calls
% * check_minimality (embedded)
+% * minrealold (embedded)
% =========================================================================
-% Copyright (C) 2019 Dynare Team
+% Copyright (C) 2019-2020 Dynare Team
%
% This file is part of Dynare.
%
@@ -60,102 +76,123 @@ function [CheckCO,minnx,minA,minB,minC,minD,dminA,dminB,dminC,dminD] = get_minim
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
% =========================================================================
-
-nx = size(A,2);
-ny = size(C,1);
-nu = size(B,2);
+if nargin == 1
+ derivs_flag = 0;
+end
+realsmall = 1e-7;
+[nspred,exo_nbr] = size(SYS.B);
+varobs_nbr = size(SYS.C,1);
% Check controllability and observability conditions for full state vector
-CheckCO = check_minimality(A,B,C);
-
-if CheckCO == 1 % If model is already minimal
- minnx = nx;
- minA = A;
- minB = B;
- minC = C;
- minD = D;
- if nargout > 6
- dminA = dA;
- dminB = dB;
- dminC = dC;
- dminD = dD;
- end
+CheckCO = check_minimality(SYS.A,SYS.B,SYS.C);
+if CheckCO == 1 % If model is already minimal, we are finished
+ minns = nspred;
+ minSYS = SYS;
else
- %Model is not minimal
- realsmall = 1e-7;
- % create indices for unnecessary states
- endogstateindex = find(abs(sum(A,1))<realsmall);
- exogstateindex = find(abs(sum(A,1))>realsmall);
- minnx = length(exogstateindex);
- % remove unnecessary states from solution matrices
- A_1 = A(endogstateindex,exogstateindex);
- A_2 = A(exogstateindex,exogstateindex);
- B_2 = B(exogstateindex,:);
- C_1 = C(:,endogstateindex);
- C_2 = C(:,exogstateindex);
- D = D;
- ind_A1 = endogstateindex;
- ind_A2 = exogstateindex;
- % minimal representation
- minA = A_2;
- minB = B_2;
- minC = C_2;
- minD = D;
- % Check controllability and observability conditions
- CheckCO = check_minimality(minA,minB,minC);
-
- if CheckCO ~=1
- j=1;
- while (CheckCO==0 && j<minnx)
- combis = nchoosek(1:minnx,j);
- i=1;
- while i <= size(combis,1)
- ind_A2 = exogstateindex;
- ind_A1 = [endogstateindex ind_A2(combis(j,:))];
- ind_A2(combis(j,:)) = [];
- % remove unnecessary states from solution matrices
- A_1 = A(ind_A1,ind_A2);
- A_2 = A(ind_A2,ind_A2);
- B_2 = B(ind_A2,:);
- C_1 = C(:,ind_A1);
- C_2 = C(:,ind_A2);
- D = D;
- % minimal representation
- minA = A_2;
- minB = B_2;
- minC = C_2;
- minD = D;
- % Check controllability and observability conditions
- CheckCO = check_minimality(minA,minB,minC);
- if CheckCO == 1
- minnx = length(ind_A2);
- break;
+ %Model is not minimal
+ try
+ minreal_flag = 1;
+ % In future we will use SLICOT TB01PD.f mex file [to do @wmutschl], currently use workaround
+ [mSYS,U] = minrealold(SYS,realsmall);
+ minns = size(mSYS.A,1);
+ CheckCO = check_minimality(mSYS.A,mSYS.B,mSYS.C);
+ if CheckCO
+ minSYS.A = mSYS.A;
+ minSYS.B = mSYS.B;
+ minSYS.C = mSYS.C;
+ minSYS.D = mSYS.D;
+ if derivs_flag
+ totparam_nbr = size(SYS.dA,3);
+ minSYS.dA = zeros(minns,minns,totparam_nbr);
+ minSYS.dB = zeros(minns,exo_nbr,totparam_nbr);
+ minSYS.dC = zeros(varobs_nbr,minns,totparam_nbr);
+ % Note that orthogonal matrix U is such that (U*dA*U',U*dB,dC*U') is a Kalman decomposition of (dA,dB,dC) %
+ for jp=1:totparam_nbr
+ dA_tmp = U*SYS.dA(:,:,jp)*U';
+ dB_tmp = U*SYS.dB(:,:,jp);
+ dC_tmp = SYS.dC(:,:,jp)*U';
+ minSYS.dA(:,:,jp) = dA_tmp(1:minns,1:minns);
+ minSYS.dB(:,:,jp) = dB_tmp(1:minns,:);
+ minSYS.dC(:,:,jp) = dC_tmp(:,1:minns);
end
- i=i+1;
+ minSYS.dD = SYS.dD;
end
- j=j+1;
end
+ catch
+ minreal_flag = 0; % if something went wrong use below procedure
end
- if nargout > 6
- dminA = dA(ind_A2,ind_A2,:);
- dminB = dB(ind_A2,:,:);
- dminC = dC(:,ind_A2,:);
- dminD = dD;
+
+ if minreal_flag == 0
+ fprintf('Use naive brute-force approach to find minimal state space system.\n These computations may be inaccurate/wrong as ''minreal'' is not available, please raise an issue on GitLab or the forum\n')
+ % create indices for unnecessary states
+ exogstateindex = find(abs(sum(SYS.A,1))>realsmall);
+ minns = length(exogstateindex);
+ % remove unnecessary states from solution matrices
+ A_2 = SYS.A(exogstateindex,exogstateindex);
+ B_2 = SYS.B(exogstateindex,:);
+ C_2 = SYS.C(:,exogstateindex);
+ D = SYS.D;
+ ind_A2 = exogstateindex;
+ % minimal representation
+ minSYS.A = A_2;
+ minSYS.B = B_2;
+ minSYS.C = C_2;
+ minSYS.D = D;
+
+ % Check controllability and observability conditions
+ CheckCO = check_minimality(minSYS.A,minSYS.B,minSYS.C);
+
+ if CheckCO ~=1
+ % do brute-force search
+ j=1;
+ while (CheckCO==0 && j<minns)
+ combis = nchoosek(1:minns,j);
+ i=1;
+ while i <= size(combis,1)
+ ind_A2 = exogstateindex;
+ ind_A2(combis(j,:)) = [];
+ % remove unnecessary states from solution matrices
+ A_2 = SYS.A(ind_A2,ind_A2);
+ B_2 = SYS.B(ind_A2,:);
+ C_2 = SYS.C(:,ind_A2);
+ D = SYS.D;
+ % minimal representation
+ minSYS.A = A_2;
+ minSYS.B = B_2;
+ minSYS.C = C_2;
+ minSYS.D = D;
+ % Check controllability and observability conditions
+ CheckCO = check_minimality(minSYS.A,minSYS.B,minSYS.C);
+ if CheckCO == 1
+ minns = length(ind_A2);
+ break;
+ end
+ i=i+1;
+ end
+ j=j+1;
+ end
+ end
+ if derivs_flag
+ minSYS.dA = SYS.dA(ind_A2,ind_A2,:);
+ minSYS.dB = SYS.dB(ind_A2,:,:);
+ minSYS.dC = SYS.dC(:,ind_A2,:);
+ minSYS.dD = SYS.dD;
+ end
end
end
-function CheckCO = check_minimality(A,B,C)
+function CheckCO = check_minimality(a,b,c)
%% This function computes the controllability and the observability matrices of the ABCD system and checks if the system is minimal
%
% Inputs: Solution matrices A,B,C of ABCD representation of a DSGE model
% Outputs: CheckCO: equals 1, if both rank conditions for observability and controllability are fullfilled
-n = size(A,1);
-CC = B; % Initialize controllability matrix
-OO = C; % Initialize observability matrix
+n = size(a,1);
+CC = b; % Initialize controllability matrix
+OO = c; % Initialize observability matrix
if n >= 2
for indn = 1:1:n-1
- CC = [CC, (A^indn)*B]; % Set up controllability matrix
- OO = [OO; C*(A^indn)]; % Set up observability matrix
+ CC = [CC, (a^indn)*b]; % Set up controllability matrix
+ OO = [OO; c*(a^indn)]; % Set up observability matrix
end
end
CheckC = (rank(full(CC))==n); % Check rank of controllability matrix
@@ -163,4 +200,87 @@ CheckO = (rank(full(OO))==n); % Check rank of observability matrix
CheckCO = CheckC&CheckO; % equals 1 if minimal state
end % check_minimality end
+function [mSYS,U] = minrealold(SYS,tol)
+ % This is a temporary replacement for minreal, will be replaced by a mex file from SLICOT TB01PD.f soon
+ a = SYS.A;
+ b = SYS.B;
+ c = SYS.C;
+ [ns,nu] = size(b);
+ [am,bm,cm,U,k] = ControllabilityStaircaseRosenbrock(a,b,c,tol);
+ kk = sum(k);
+ nu = ns - kk;
+ nn = nu;
+ am = am(nu+1:ns,nu+1:ns);
+ bm = bm(nu+1:ns,:);
+ cm = cm(:,nu+1:ns);
+ ns = ns - nu;
+ if ns
+ [am,bm,cm,t,k] = ObservabilityStaircaseRosenbrock(am,bm,cm,tol);
+ kk = sum(k);
+ nu = ns - kk;
+ nn = nn + nu;
+ am = am(nu+1:ns,nu+1:ns);
+ bm = bm(nu+1:ns,:);
+ cm = cm(:,nu+1:ns);
+ end
+ mSYS.A = am;
+ mSYS.B = bm;
+ mSYS.C = cm;
+ mSYS.D = SYS.D;
+end
+
+
+function [abar,bbar,cbar,t,k] = ObservabilityStaircaseRosenbrock(a,b,c,tol)
+ %Observability staircase form
+ [aa,bb,cc,t,k] = ControllabilityStaircaseRosenbrock(a',c',b',tol);
+ abar = aa'; bbar = cc'; cbar = bb';
+end
+
+function [abar,bbar,cbar,t,k] = ControllabilityStaircaseRosenbrock(a, b, c, tol)
+ % Controllability staircase algorithm of Rosenbrock, 1968
+ [ra,ca] = size(a);
+ [rb,cb] = size(b);
+ ptjn1 = eye(ra);
+ ajn1 = a;
+ bjn1 = b;
+ rojn1 = cb;
+ deltajn1 = 0;
+ sigmajn1 = ra;
+ k = zeros(1,ra);
+ if nargin == 3
+ tol = ra*norm(a,1)*eps;
+ end
+ for jj = 1 : ra
+ [uj,sj,vj] = svd(bjn1);
+ [rsj,csj] = size(sj);
+ p = flip(eye(rsj),2);
+ p = permute(p,[2 1 3:ndims(eye(rsj))]);
+ uj = uj*p;
+ bb = uj'*bjn1;
+ roj = rank(bb,tol);
+ [rbb,cbb] = size(bb);
+ sigmaj = rbb - roj;
+ sigmajn1 = sigmaj;
+ k(jj) = roj;
+ if roj == 0, break, end
+ if sigmaj == 0, break, end
+ abxy = uj' * ajn1 * uj;
+ aj = abxy(1:sigmaj,1:sigmaj);
+ bj = abxy(1:sigmaj,sigmaj+1:sigmaj+roj);
+ ajn1 = aj;
+ bjn1 = bj;
+ [ruj,cuj] = size(uj);
+ ptj = ptjn1 * ...
+ [uj zeros(ruj,deltajn1); ...
+ zeros(deltajn1,cuj) eye(deltajn1)];
+ ptjn1 = ptj;
+ deltaj = deltajn1 + roj;
+ deltajn1 = deltaj;
+ end
+ t = ptjn1';
+ abar = t * a * t';
+ bbar = t * b;
+ cbar = c * t';
+end
+
end % Main function end
diff --git a/tests/minimal state space system/as2007_minimal.mod b/tests/minimal state space system/as2007_minimal.mod
new file mode 100644
index 0000000000000000000000000000000000000000..788a70ca02c132621cdc8586e14b4c7760da9ab4
--- /dev/null
+++ b/tests/minimal state space system/as2007_minimal.mod
@@ -0,0 +1,100 @@
+% this is the Smets and Wouters (2007) model for which Komunjer and Ng (2011)
+% derived the minimal state space system. In Dynare, however, we use more
+% powerful minreal function
+% created by Willi Mutschler (@wmutschl, willi@mutschler.eu)
+% =========================================================================
+% Copyright (C) 2020 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 <http://www.gnu.org/licenses/>.
+% =========================================================================
+
+var y R g z c dy p YGR INFL INT;
+varobs y R p c YGR INFL INT;
+varexo e_r e_g e_z;
+parameters tau phi psi1 psi2 rhor rhog rhoz rrst pist gamst nu cyst;
+
+rrst = 1.0000;
+pist = 3.2000;
+gamst= 0.5500;
+tau = 2.0000;
+nu = 0.1000;
+kap = 0.3300;
+phi = tau*(1-nu)/nu/kap/exp(pist/400)^2;
+cyst = 0.8500;
+psi1 = 1.5000;
+psi2 = 0.1250;
+rhor = 0.7500;
+rhog = 0.9500;
+rhoz = 0.9000;
+
+
+model;
+#pist2 = exp(pist/400);
+#rrst2 = exp(rrst/400);
+#bet = 1/rrst2;
+#gst = 1/cyst;
+#cst = (1-nu)^(1/tau);
+#yst = cst*gst;
+1 = exp(-tau*c(+1)+tau*c+R-z(+1)-p(+1));
+(1-nu)/nu/phi/(pist2^2)*(exp(tau*c)-1) = (exp(p)-1)*((1-1/2/nu)*exp(p)+1/2/nu) - bet*(exp(p(+1))-1)*exp(-tau*c(+1)+tau*c+dy(+1)+p(+1));
+exp(c-y) = exp(-g) - phi*pist2^2*gst/2*(exp(p)-1)^2;
+R = rhor*R(-1) + (1-rhor)*psi1*p + (1-rhor)*psi2*(y-g) + e_r;
+g = rhog*g(-1) + e_g;
+z = rhoz*z(-1) + e_z;
+YGR = gamst+100*(dy+z);
+INFL = pist+400*p;
+INT = pist+rrst+4*gamst+400*R;
+dy = y - y(-1);
+end;
+
+shocks;
+var e_r; stderr 0.2/100;
+var e_g; stderr 0.6/100;
+var e_z; stderr 0.3/100;
+end;
+
+steady_state_model;
+z=0; g=0; c=0; y=0; p=0; R=0; dy=0;
+YGR=gamst; INFL=pist; INT=pist+rrst+4*gamst;
+end;
+stoch_simul(order=1,irf=0,periods=0);
+options_.qz_criterium = 1;
+
+indx = [M_.nstatic+(1:M_.nspred)]';
+indy = 1:M_.endo_nbr';
+
+SS.A = oo_.dr.ghx(indx,:);
+SS.B = oo_.dr.ghu(indx,:);
+SS.C = oo_.dr.ghx(indy,:);
+SS.D = oo_.dr.ghu(indy,:);
+
+[CheckCO,minnx,minSS] = get_minimal_state_representation(SS,0);
+
+Sigmax_full = lyapunov_symm(SS.A, SS.B*M_.Sigma_e*SS.B', options_.lyapunov_fixed_point_tol, options_.qz_criterium, options_.lyapunov_complex_threshold, 1, options_.debug);
+Sigmay_full = SS.C*Sigmax_full*SS.C' + SS.D*M_.Sigma_e*SS.D';
+
+Sigmax_min = lyapunov_symm(minSS.A, minSS.B*M_.Sigma_e*minSS.B', options_.lyapunov_fixed_point_tol, options_.qz_criterium, options_.lyapunov_complex_threshold, 1, options_.debug);
+Sigmay_min = minSS.C*Sigmax_min*minSS.C' + minSS.D*M_.Sigma_e*minSS.D';
+
+([Sigmay_full(:) - Sigmay_min(:)]')
+sqrt(([diag(Sigmay_full), diag(Sigmay_min)]'))
+dx = norm( Sigmay_full - Sigmay_min, Inf);
+if dx > 1e-12
+ error('something wrong with minimal state space computations')
+else
+ fprintf('numerical error for moments computed from minimal state system is %d\n',dx)
+end
+
diff --git a/tests/minimal state space system/sw_minimal.mod b/tests/minimal state space system/sw_minimal.mod
new file mode 100644
index 0000000000000000000000000000000000000000..2ced78e1a0d09f1bb2830745b5aecf89a299136e
--- /dev/null
+++ b/tests/minimal state space system/sw_minimal.mod
@@ -0,0 +1,432 @@
+% this is the Smets and Wouters (2007) model for which Komunjer and Ng (2011)
+% derived the minimal state space system. In Dynare, however, we use more
+% powerful minreal function
+/*
+ * This file provides replication files for
+ * Smets, Frank and Wouters, Rafael (2007): "Shocks and Frictions in US Business Cycles: A Bayesian
+ * DSGE Approach", American Economic Review, 97(3), 586-606, that are compatible with Dynare 4.5 onwards
+ *
+ * To replicate the full results, you have to get back to the original replication files available at
+ * https://www.aeaweb.org/articles.php?doi=10.1257/aer.97.3.586 and include the respective estimation commands and mode-files.
+ *
+ * Notes:
+ * - The consumption Euler equation in the paper, equation (2), premultiplies the risk premium process \varepsilon_t^b,
+ * denoted by b in this code, by the coefficient c_3. In the code this prefactor is omitted by setting the
+ * coefficient to 1. As a consequence, b in this code actually is b:=c_3*\varepsilon_t^b. As a consequence, in
+ * the arbitrage equation for the value of capital in the paper, equation (4), the term 1*\varepsilon_t^b
+ * is replaced by 1/c_3*b, which is equal to \varepsilon_t^b given the above redefinition. This rescaling also explains why the
+ * standard deviation of the risk premium shock in the AR(1)-process for b has a different standard deviation than reported
+ * in the paper. However, the results are unaffected by this scaling factor (except for the fact that the posterior distribution
+ * reported in the paper cannot be directly translated to the present mod-file due to parameter correlation in the posterior.
+ * - As pointed out in Del Negro/Schorfheide (2012): "Notes on New-Keynesian Models"
+ * in the code implementation of equation (8) for both the flex price and the sticky price/wage economy,
+ * there is a (1/(1+cbetabar*cgamma)) missing in the i_2 in front of q_t (denoted qs in the code).
+ * Equation (8) in the paper reads:
+ * (1-(1-delta)/gamma)*(1+beta*gamma^(1-sigma))*gamma^2*varphi
+ * which translates to the code snippet:
+ * (1-(1-ctou)/cgamma)*(1+cbetabar*cgamma)*cgamma^2*csadjcost
+ * But the code implements
+ * (1-(1-ctou)/cgamma)*cgamma^2*csadjcost
+ * which corresponds to an equation reading
+ * (1-(1-delta)/gamma)*gamma^2*varphi
+ * - Chib/Ramamurthy (2010): "Tailored randomized block MCMC methods with application to DSGE models", Journal of Econometrics, 155, pp. 19-38
+ * have pointed out that the mode reported in the original Smets/Wouters (2007) paper is not actually the mode. \bar \pi (constepinf) is estimated lower
+ * while \bar \l (constelab) is higher.
+ * - Note that at the prior mean, [cmap,crhopinf] and [cmaw,crhow] are pairwise collinear. Thus, running identification at the prior
+ * mean will return a warning. But this is only a local issue. These parameters are only indistinguishable at the prior mean, but not
+ * at different points.
+ * - In the prior Table 1A in the paper, the
+ * - habit parameter $\lambda$ is erroneously labeled h
+ * - the fixed cost parameter $\phi_p$ is labeled $\Phi$
+ * - Table 1B claims that $\rho_{ga}$ follows a beta prior with B(0.5,0.2^2), but the code shows that it actually
+ * follows a normal distribution with N(0.5,0.25^2)
+ *
+ * This file was originally written by Frank Smets and Rafeal Wouters and has been updated by
+ * Johannes Pfeifer.
+ *
+ * Please note that the following copyright notice only applies to this Dynare
+ * implementation of the model
+ */
+
+/*
+ * Copyright (C) 2007-2013 Frank Smets and Raf Wouters
+ * Copyright (C) 2013-15 Johannes Pfeifer
+ * Copyright (C) 2020 Dynare Team
+ *
+ * This 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.
+ *
+ * This file 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 can receive a copy of the GNU General Public License
+ * at <http://www.gnu.org/licenses/>.
+ */
+
+var labobs ${lHOURS}$ (long_name='log hours worked')
+ robs ${FEDFUNDS}$ (long_name='Federal funds rate')
+ pinfobs ${dlP}$ (long_name='Inflation')
+ dy ${dlGDP}$ (long_name='Output growth rate')
+ dc ${dlCONS}$ (long_name='Consumption growth rate')
+ dinve ${dlINV}$ (long_name='Investment growth rate')
+ dw ${dlWAG}$ (long_name='Wage growth rate')
+ ewma ${\eta^{w,aux}}$ (long_name='Auxiliary wage markup moving average variable')
+ epinfma ${\eta^{p,aux}}$ (long_name='Auxiliary price markup moving average variable')
+ zcapf ${z^{flex}}$ (long_name='Capital utilization rate flex price economy')
+ rkf ${r^{k,flex}}$ (long_name='rental rate of capital flex price economy')
+ kf ${k^{s,flex}}$ (long_name='Capital services flex price economy')
+ pkf ${q^{flex}}$ (long_name='real value of existing capital stock flex price economy')
+ cf ${c^{flex}}$ (long_name='Consumption flex price economy')
+ invef ${i^{flex}}$ (long_name='Investment flex price economy')
+ yf ${y^{flex}}$ (long_name='Output flex price economy')
+ labf ${l^{flex}}$ (long_name='hours worked flex price economy')
+ wf ${w^{flex}}$ (long_name='real wage flex price economy')
+ rrf ${r^{flex}}$ (long_name='real interest rate flex price economy')
+ mc ${\mu_p}$ (long_name='gross price markup')
+ zcap ${z}$ (long_name='Capital utilization rate')
+ rk ${r^{k}}$ (long_name='rental rate of capital')
+ k ${k^{s}}$ (long_name='Capital services')
+ pk ${q}$ (long_name='real value of existing capital stock')
+ c ${c}$ (long_name='Consumption')
+ inve ${i}$ (long_name='Investment')
+ y ${y}$ (long_name='Output')
+ lab ${l}$ (long_name='hours worked')
+ pinf ${\pi}$ (long_name='Inflation')
+ w ${w}$ (long_name='real wage')
+ r ${r}$ (long_name='nominal interest rate')
+ a ${\varepsilon_a}$ (long_name='productivity process')
+ b ${c_2*\varepsilon_t^b}$ (long_name='Scaled risk premium shock')
+ g ${\varepsilon^g}$ (long_name='Exogenous spending')
+ qs ${\varepsilon^i}$ (long_name='Investment-specific technology')
+ ms ${\varepsilon^r}$ (long_name='Monetary policy shock process')
+ spinf ${\varepsilon^p}$ (long_name='Price markup shock process')
+ sw ${\varepsilon^w}$ (long_name='Wage markup shock process')
+ kpf ${k^{flex}}$ (long_name='Capital stock flex price economy')
+ kp ${k}$ (long_name='Capital stock')
+ ;
+
+varexo ea ${\eta^a}$ (long_name='productivity shock')
+ eb ${\eta^b}$ (long_name='Investment-specific technology shock')
+ eg ${\eta^g}$ (long_name='Spending shock')
+ eqs ${\eta^i}$ (long_name='Investment-specific technology shock')
+ em ${\eta^m}$ (long_name='Monetary policy shock')
+ epinf ${\eta^{p}}$ (long_name='Price markup shock')
+ ew ${\eta^{w}}$ (long_name='Wage markup shock')
+ ;
+
+parameters curvw ${\varepsilon_w}$ (long_name='Curvature Kimball aggregator wages')
+ cgy ${\rho_{ga}}$ (long_name='Feedback technology on exogenous spending')
+ curvp ${\varepsilon_p}$ (long_name='Curvature Kimball aggregator prices')
+ constelab ${\bar l}$ (long_name='steady state hours')
+ constepinf ${\bar \pi}$ (long_name='steady state inflation rate')
+ constebeta ${100(\beta^{-1}-1)}$ (long_name='time preference rate in percent')
+ cmaw ${\mu_w}$ (long_name='coefficient on MA term wage markup')
+ cmap ${\mu_p}$ (long_name='coefficient on MA term price markup')
+ calfa ${\alpha}$ (long_name='capital share')
+ czcap ${\psi}$ (long_name='capacity utilization cost')
+ csadjcost ${\varphi}$ (long_name='investment adjustment cost')
+ ctou ${\delta}$ (long_name='depreciation rate')
+ csigma ${\sigma_c}$ (long_name='risk aversion')
+ chabb ${\lambda}$ (long_name='external habit degree')
+ ccs ${d_4}$ (long_name='Unused parameter')
+ cinvs ${d_3}$ (long_name='Unused parameter')
+ cfc ${\phi_p}$ (long_name='fixed cost share')
+ cindw ${\iota_w}$ (long_name='Indexation to past wages')
+ cprobw ${\xi_w}$ (long_name='Calvo parameter wages')
+ cindp ${\iota_p}$ (long_name='Indexation to past prices')
+ cprobp ${\xi_p}$ (long_name='Calvo parameter prices')
+ csigl ${\sigma_l}$ (long_name='Frisch elasticity')
+ clandaw ${\phi_w}$ (long_name='Gross markup wages')
+ crdpi ${r_{\Delta \pi}}$ (long_name='Unused parameter')
+ crpi ${r_{\pi}}$ (long_name='Taylor rule inflation feedback')
+ crdy ${r_{\Delta y}}$ (long_name='Taylor rule output growth feedback')
+ cry ${r_{y}}$ (long_name='Taylor rule output level feedback')
+ crr ${\rho}$ (long_name='interest rate persistence')
+ crhoa ${\rho_a}$ (long_name='persistence productivity shock')
+ crhoas ${d_2}$ (long_name='Unused parameter')
+ crhob ${\rho_b}$ (long_name='persistence risk premium shock')
+ crhog ${\rho_g}$ (long_name='persistence spending shock')
+ crhols ${d_1}$ (long_name='Unused parameter')
+ crhoqs ${\rho_i}$ (long_name='persistence risk premium shock')
+ crhoms ${\rho_r}$ (long_name='persistence monetary policy shock')
+ crhopinf ${\rho_p}$ (long_name='persistence price markup shock')
+ crhow ${\rho_w}$ (long_name='persistence wage markup shock')
+ ctrend ${\bar \gamma}$ (long_name='net growth rate in percent')
+ cg ${\frac{\bar g}{\bar y}}$ (long_name='steady state exogenous spending share')
+ ;
+
+// fixed parameters
+ctou=.025;
+clandaw=1.5;
+cg=0.18;
+curvp=10;
+curvw=10;
+
+// estimated parameters initialisation
+calfa=.24;
+cbeta=.9995;
+csigma=1.5;
+cfc=1.5;
+cgy=0.51;
+
+csadjcost= 6.0144;
+chabb= 0.6361;
+cprobw= 0.8087;
+csigl= 1.9423;
+cprobp= 0.6;
+cindw= 0.3243;
+cindp= 0.47;
+czcap= 0.2696;
+crpi= 1.488;
+crr= 0.8762;
+cry= 0.0593;
+crdy= 0.2347;
+
+crhoa= 0.9977;
+crhob= 0.5799;
+crhog= 0.9957;
+crhols= 0.9928;
+crhoqs= 0.7165;
+crhoas=1;
+crhoms=0;
+crhopinf=0;
+crhow=0;
+cmap = 0;
+cmaw = 0;
+
+constelab=0;
+
+%% Added by JP to provide full calibration of model before estimation
+constepinf=0.7;
+constebeta=0.7420;
+ctrend=0.3982;
+
+model(linear);
+//deal with parameter dependencies; taken from usmodel_stst.mod
+#cpie=1+constepinf/100; %gross inflation rate
+#cgamma=1+ctrend/100 ; %gross growth rate
+#cbeta=1/(1+constebeta/100); %discount factor
+
+#clandap=cfc; %fixed cost share/gross price markup
+#cbetabar=cbeta*cgamma^(-csigma); %growth-adjusted discount factor in Euler equation
+#cr=cpie/(cbeta*cgamma^(-csigma)); %steady state net real interest rate
+#crk=(cbeta^(-1))*(cgamma^csigma) - (1-ctou); %R^k_{*}: steady state rental rate
+#cw = (calfa^calfa*(1-calfa)^(1-calfa)/(clandap*crk^calfa))^(1/(1-calfa)); %steady state real wage
+//cw = (calfa^calfa*(1-calfa)^(1-calfa)/(clandap*((cbeta^(-1))*(cgamma^csigma) - (1-ctou))^calfa))^(1/(1-calfa));
+#cikbar=(1-(1-ctou)/cgamma); %k_1 in equation LOM capital, equation (8)
+#cik=(1-(1-ctou)/cgamma)*cgamma; %i_k: investment-capital ratio
+#clk=((1-calfa)/calfa)*(crk/cw); %labor to capital ratio
+#cky=cfc*(clk)^(calfa-1); %k_y: steady state output ratio
+#ciy=cik*cky; %consumption-investment ratio
+#ccy=1-cg-cik*cky; %consumption-output ratio
+#crkky=crk*cky; %z_y=R_{*}^k*k_y
+#cwhlc=(1/clandaw)*(1-calfa)/calfa*crk*cky/ccy; %W^{h}_{*}*L_{*}/C_{*} used in c_2 in equation (2)
+#cwly=1-crk*cky; %unused parameter
+
+#conster=(cr-1)*100; %steady state federal funds rate ($\bar r$)
+
+// flexible economy
+
+ [name='FOC labor with mpl expressed as function of rk and w, flex price economy']
+ 0*(1-calfa)*a + 1*a = calfa*rkf+(1-calfa)*(wf) ;
+ [name='FOC capacity utilization, flex price economy']
+ zcapf = (1/(czcap/(1-czcap)))* rkf ;
+ [name='Firm FOC capital, flex price economy']
+ rkf = (wf)+labf-kf ;
+ [name='Definition capital services, flex price economy']
+ kf = kpf(-1)+zcapf ;
+ [name='Investment Euler Equation, flex price economy']
+ invef = (1/(1+cbetabar*cgamma))* ( invef(-1) + cbetabar*cgamma*invef(1)+(1/(cgamma^2*csadjcost))*pkf ) +qs ;
+ [name='Arbitrage equation value of capital, flex price economy']
+ pkf = -rrf-0*b+(1/((1-chabb/cgamma)/(csigma*(1+chabb/cgamma))))*b +(crk/(crk+(1-ctou)))*rkf(1) + ((1-ctou)/(crk+(1-ctou)))*pkf(1) ;
+ [name='Consumption Euler Equation, flex price economy']
+ cf = (chabb/cgamma)/(1+chabb/cgamma)*cf(-1) + (1/(1+chabb/cgamma))*cf(+1) +((csigma-1)*cwhlc/(csigma*(1+chabb/cgamma)))*(labf-labf(+1)) - (1-chabb/cgamma)/(csigma*(1+chabb/cgamma))*(rrf+0*b) + b ;
+ [name='Aggregate Resource Constraint, flex price economy']
+ yf = ccy*cf+ciy*invef+g + crkky*zcapf ;
+ [name='Aggregate Production Function, flex price economy']
+ yf = cfc*( calfa*kf+(1-calfa)*labf +a );
+ [name='Wage equation, flex price economy']
+ wf = csigl*labf +(1/(1-chabb/cgamma))*cf - (chabb/cgamma)/(1-chabb/cgamma)*cf(-1) ;
+ [name='Law of motion for capital, flex price economy (see header notes)']
+ kpf = (1-cikbar)*kpf(-1)+(cikbar)*invef + (cikbar)*(cgamma^2*csadjcost)*qs ;
+
+// sticky price - wage economy
+ [name='FOC labor with mpl expressed as function of rk and w, SW Equation (9)']
+ mc = calfa*rk+(1-calfa)*(w) - 1*a - 0*(1-calfa)*a ;
+ [name='FOC capacity utilization, SW Equation (7)']
+ zcap = (1/(czcap/(1-czcap)))* rk ;
+ [name='Firm FOC capital, SW Equation (11)']
+ rk = w+lab-k ;
+ [name='Definition capital services, SW Equation (6)']
+ k = kp(-1)+zcap ;
+ [name='Investment Euler Equation, SW Equation (3)']
+ inve = (1/(1+cbetabar*cgamma))* (inve(-1) + cbetabar*cgamma*inve(1)+(1/(cgamma^2*csadjcost))*pk ) +qs ;
+ [name='Arbitrage equation value of capital, SW Equation (4)']
+ pk = -r+pinf(1)-0*b
+ + (1/((1-chabb/cgamma)/(csigma*(1+chabb/cgamma))))*b
+ + (crk/(crk+(1-ctou)))*rk(1)
+ + ((1-ctou)/(crk+(1-ctou)))*pk(1) ;
+ [name='Consumption Euler Equation, SW Equation (2)']
+ c = (chabb/cgamma)/(1+chabb/cgamma)*c(-1)
+ + (1/(1+chabb/cgamma))*c(+1)
+ +((csigma-1)*cwhlc/(csigma*(1+chabb/cgamma)))*(lab-lab(+1))
+ - (1-chabb/cgamma)/(csigma*(1+chabb/cgamma))*(r-pinf(+1) + 0*b) +b ;
+ [name='Aggregate Resource Constraint, SW Equation (1)']
+ y = ccy*c+ciy*inve+g + 1*crkky*zcap ;
+ [name='Aggregate Production Function, SW Equation (5)']
+ y = cfc*( calfa*k+(1-calfa)*lab +a );
+ [name='New Keynesian Phillips Curve, SW Equation (10)']
+ pinf = (1/(1+cbetabar*cgamma*cindp)) * ( cbetabar*cgamma*pinf(1) +cindp*pinf(-1)
+ +((1-cprobp)*(1-cbetabar*cgamma*cprobp)/cprobp)/((cfc-1)*curvp+1)*(mc) ) + spinf ;
+ [name='Wage Phillips Curve, SW Equation (13), with (12) plugged for mu_w']
+ w = (1/(1+cbetabar*cgamma))*w(-1)
+ +(cbetabar*cgamma/(1+cbetabar*cgamma))*w(1)
+ +(cindw/(1+cbetabar*cgamma))*pinf(-1)
+ -(1+cbetabar*cgamma*cindw)/(1+cbetabar*cgamma)*pinf
+ +(cbetabar*cgamma)/(1+cbetabar*cgamma)*pinf(1)
+ +(1-cprobw)*(1-cbetabar*cgamma*cprobw)/((1+cbetabar*cgamma)*cprobw)*(1/((clandaw-1)*curvw+1))*
+ (csigl*lab + (1/(1-chabb/cgamma))*c - ((chabb/cgamma)/(1-chabb/cgamma))*c(-1) -w)
+ + 1*sw ;
+ [name='Taylor rule, SW Equation (14)']
+ r = crpi*(1-crr)*pinf
+ +cry*(1-crr)*(y-yf)
+ +crdy*(y-yf-y(-1)+yf(-1))
+ +crr*r(-1)
+ +ms ;
+ [name='Law of motion for productivity']
+ a = crhoa*a(-1) + ea;
+ [name='Law of motion for risk premium']
+ b = crhob*b(-1) + eb;
+ [name='Law of motion for spending process']
+ g = crhog*(g(-1)) + eg + cgy*ea;
+ [name='Law of motion for investment specific technology shock process']
+ qs = crhoqs*qs(-1) + eqs;
+ [name='Law of motion for monetary policy shock process']
+ ms = crhoms*ms(-1) + em;
+ [name='Law of motion for price markup shock process']
+ spinf = crhopinf*spinf(-1) + epinfma - cmap*epinfma(-1);
+ epinfma=epinf;
+ [name='Law of motion for wage markup shock process']
+ sw = crhow*sw(-1) + ewma - cmaw*ewma(-1) ;
+ ewma=ew;
+ [name='Law of motion for capital, SW Equation (8) (see header notes)']
+ kp = (1-cikbar)*kp(-1)+cikbar*inve + cikbar*cgamma^2*csadjcost*qs ;
+
+// measurement equations
+[name='Observation equation output']
+dy=y-y(-1)+ctrend;
+[name='Observation equation consumption']
+dc=c-c(-1)+ctrend;
+[name='Observation equation investment']
+dinve=inve-inve(-1)+ctrend;
+[name='Observation equation real wage']
+dw=w-w(-1)+ctrend;
+[name='Observation equation inflation']
+pinfobs = 1*(pinf) + constepinf;
+[name='Observation equation interest rate']
+robs = 1*(r) + conster;
+[name='Observation equation hours worked']
+labobs = lab + constelab;
+
+end;
+
+steady_state_model;
+dy=ctrend;
+dc=ctrend;
+dinve=ctrend;
+dw=ctrend;
+pinfobs = constepinf;
+robs = (((1+constepinf/100)/((1/(1+constebeta/100))*(1+ctrend/100)^(-csigma)))-1)*100;
+labobs = constelab;
+end;
+
+shocks;
+var ea;
+stderr 0.4618;
+var eb;
+stderr 1.8513;
+var eg;
+stderr 0.6090;
+var eqs;
+stderr 0.6017;
+var em;
+stderr 0.2397;
+var epinf;
+stderr 0.1455;
+var ew;
+stderr 0.2089;
+end;
+
+
+estimated_params;
+// PARAM NAME, INITVAL, LB, UB, PRIOR_SHAPE, PRIOR_P1, PRIOR_P2, PRIOR_P3, PRIOR_P4, JSCALE
+// PRIOR_SHAPE: BETA_PDF, GAMMA_PDF, NORMAL_PDF, INV_GAMMA_PDF
+stderr ea,0.4618,0.01,3,INV_GAMMA_PDF,0.1,2;
+stderr eb,0.1818513,0.025,5,INV_GAMMA_PDF,0.1,2;
+stderr eg,0.6090,0.01,3,INV_GAMMA_PDF,0.1,2;
+stderr eqs,0.46017,0.01,3,INV_GAMMA_PDF,0.1,2;
+stderr em,0.2397,0.01,3,INV_GAMMA_PDF,0.1,2;
+stderr epinf,0.1455,0.01,3,INV_GAMMA_PDF,0.1,2;
+stderr ew,0.2089,0.01,3,INV_GAMMA_PDF,0.1,2;
+crhoa,.9676 ,.01,.9999,BETA_PDF,0.5,0.20;
+crhob,.2703,.01,.9999,BETA_PDF,0.5,0.20;
+crhog,.9930,.01,.9999,BETA_PDF,0.5,0.20;
+crhoqs,.5724,.01,.9999,BETA_PDF,0.5,0.20;
+crhoms,.3,.01,.9999,BETA_PDF,0.5,0.20;
+crhopinf,.8692,.01,.9999,BETA_PDF,0.5,0.20;
+crhow,.9546,.001,.9999,BETA_PDF,0.5,0.20;
+cmap,.7652,0.01,.9999,BETA_PDF,0.5,0.2;
+cmaw,.8936,0.01,.9999,BETA_PDF,0.5,0.2;
+csadjcost,6.3325,2,15,NORMAL_PDF,4,1.5;
+csigma,1.2312,0.25,3,NORMAL_PDF,1.50,0.375;
+chabb,0.7205,0.001,0.99,BETA_PDF,0.7,0.1;
+cprobw,0.7937,0.3,0.95,BETA_PDF,0.5,0.1;
+csigl,2.8401,0.25,10,NORMAL_PDF,2,0.75;
+cprobp,0.7813,0.5,0.95,BETA_PDF,0.5,0.10;
+cindw,0.4425,0.01,0.99,BETA_PDF,0.5,0.15;
+cindp,0.3291,0.01,0.99,BETA_PDF,0.5,0.15;
+czcap,0.2648,0.01,1,BETA_PDF,0.5,0.15;
+cfc,1.4672,1.0,3,NORMAL_PDF,1.25,0.125;
+crpi,1.7985,1.0,3,NORMAL_PDF,1.5,0.25;
+crr,0.8258,0.5,0.975,BETA_PDF,0.75,0.10;
+cry,0.0893,0.001,0.5,NORMAL_PDF,0.125,0.05;
+crdy,0.2239,0.001,0.5,NORMAL_PDF,0.125,0.05;
+constepinf,0.7,0.1,2.0,GAMMA_PDF,0.625,0.1;//20;
+constebeta,0.7420,0.01,2.0,GAMMA_PDF,0.25,0.1;//0.20;
+constelab,1.2918,-10.0,10.0,NORMAL_PDF,0.0,2.0;
+ctrend,0.3982,0.1,0.8,NORMAL_PDF,0.4,0.10;
+cgy,0.05,0.01,2.0,NORMAL_PDF,0.5,0.25;
+calfa,0.24,0.01,1.0,NORMAL_PDF,0.3,0.05;
+end;
+
+stoch_simul(order=1,irf=0,periods=0);
+options_.qz_criterium = 1;
+
+indx = [M_.nstatic+(1:M_.nspred)]';
+indy = 1:M_.endo_nbr';
+
+SS.A = oo_.dr.ghx(indx,:);
+SS.B = oo_.dr.ghu(indx,:);
+SS.C = oo_.dr.ghx(indy,:);
+SS.D = oo_.dr.ghu(indy,:);
+
+[CheckCO,minnx,minSS] = get_minimal_state_representation(SS,0);
+
+Sigmax_full = lyapunov_symm(SS.A, SS.B*M_.Sigma_e*SS.B', options_.lyapunov_fixed_point_tol, options_.qz_criterium, options_.lyapunov_complex_threshold, 1, options_.debug);
+Sigmay_full = SS.C*Sigmax_full*SS.C' + SS.D*M_.Sigma_e*SS.D';
+
+Sigmax_min = lyapunov_symm(minSS.A, minSS.B*M_.Sigma_e*minSS.B', options_.lyapunov_fixed_point_tol, options_.qz_criterium, options_.lyapunov_complex_threshold, 1, options_.debug);
+Sigmay_min = minSS.C*Sigmax_min*minSS.C' + minSS.D*M_.Sigma_e*minSS.D';
+
+([Sigmay_full(:) - Sigmay_min(:)]')
+sqrt(([diag(Sigmay_full), diag(Sigmay_min)]'))
+dx = norm( Sigmay_full - Sigmay_min, Inf);
+if dx > 1e-8
+ error(sprintf('something wrong with minimal state space computations, as numerical error is %d',dx))
+else
+ fprintf('numerical error for moments computed from minimal state system is %d\n',dx)
+end