Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found
Select Git revision

Target

Select target project
  • giovanma/dynare
  • giorgiomas/dynare
  • Vermandel/dynare
  • Dynare/dynare
  • normann/dynare
  • MichelJuillard/dynare
  • wmutschl/dynare
  • FerhatMihoubi/dynare
  • sebastien/dynare
  • lnsongxf/dynare
  • rattoma/dynare
  • CIMERS/dynare
  • FredericKarame/dynare
  • SumuduK/dynare
  • MinjeJeon/dynare
  • camilomrch/dynare
  • DoraK/dynare
  • avtishin/dynare
  • selma/dynare
  • claudio_olguin/dynare
  • jeffjiang07/dynare
  • EthanSystem/dynare
  • stepan-a/dynare
  • wjgatt/dynare
  • JohannesPfeifer/dynare
  • gboehl/dynare
  • chskcau/dynare-doc-fixes
27 results
Select Git revision
Show changes
Commits on Source (21)
Show whitespace changes
Inline Side-by-side
Showing
with 22 additions and 1294 deletions
.gitignore
View file @ da322f81
...@@ -139,34 +139,5 @@ mex/build/matlab/run_m2html.m ...@@ -139,34 +139,5 @@ mex/build/matlab/run_m2html.m
# Windows # Windows
/windows/dynare-version.nsi /windows/dynare-version.nsi
# Estimation DLL tests
/mex/sources/estimation/tests/test-dr
/mex/sources/estimation/tests/test-dr.exe
/mex/sources/estimation/tests/testModelSolution
/mex/sources/estimation/tests/testModelSolution.exe
/mex/sources/estimation/tests/testInitKalman
/mex/sources/estimation/tests/testInitKalman.exe
/mex/sources/estimation/tests/testKalman
/mex/sources/estimation/tests/testKalman.exe
/mex/sources/estimation/tests/testPDF
/mex/sources/estimation/tests/testPDF.exe
/mex/sources/estimation/libmat/tests/test-qr
/mex/sources/estimation/libmat/tests/test-qr.exe
/mex/sources/estimation/libmat/tests/test-gsd
/mex/sources/estimation/libmat/tests/test-gsd.exe
/mex/sources/estimation/libmat/tests/test-lu
/mex/sources/estimation/libmat/tests/test-lu.exe
/mex/sources/estimation/libmat/tests/test-repmat
/mex/sources/estimation/libmat/tests/test-repmat.exe
# Misc
!/matlab/swz/c-code/Makefile
!/mex/sources/kalman/cc/Makefile
!/mex/sources/kalman/matlab/Makefile
!/mex/sources/kalman/qt/f90/Makefile
!/mex/sources/kalman/qt/test/Makefile
!/mex/sources/kalman/sylv/cc/Makefile
!/mex/sources/kalman/testing/Makefile
# MacOS stuff # MacOS stuff
.DS_Store .DS_Store
.gitmodules deleted 100644 → 0
View file @ b00b2fad
[submodule "ms-sbvar/utilities_dw"]
path = ms-sbvar/utilities_dw
url = http://www.dynare.org/git/frbatlanta/utilities_dw.git
[submodule "ms-sbvar/TZcode"]
path = ms-sbvar/TZcode
url = http://www.dynare.org/git/frbatlanta/TZcode.git
[submodule "ms-sbvar/switch_dw"]
path = ms-sbvar/switch_dw
url = http://www.dynare.org/git/frbatlanta/switch_dw.git
configure.ac
View file @ da322f81
...@@ -18,7 +18,7 @@ dnl You should have received a copy of the GNU General Public License ...@@ -18,7 +18,7 @@ dnl You should have received a copy of the GNU General Public License
dnl along with Dynare. If not, see <http://www.gnu.org/licenses/>. dnl along with Dynare. If not, see <http://www.gnu.org/licenses/>.
AC_PREREQ([2.61]) AC_PREREQ([2.61])
AC_INIT([dynare], [4.2-unstable]) AC_INIT([dynare], [4.2.0])
AC_CONFIG_SRCDIR([preprocessor/DynareMain.cc]) AC_CONFIG_SRCDIR([preprocessor/DynareMain.cc])
AM_INIT_AUTOMAKE([-Wall -Werror foreign]) AM_INIT_AUTOMAKE([-Wall -Werror foreign])
...@@ -190,10 +190,6 @@ AC_CONFIG_FILES([Makefile ...@@ -190,10 +190,6 @@ AC_CONFIG_FILES([Makefile
dynare++/kord/Makefile dynare++/kord/Makefile
dynare++/src/Makefile dynare++/src/Makefile
mex/sources/Makefile mex/sources/Makefile
mex/sources/estimation/Makefile
mex/sources/estimation/tests/Makefile
mex/sources/estimation/libmat/Makefile
mex/sources/estimation/libmat/tests/Makefile
]) ])
AC_ARG_ENABLE([matlab], AS_HELP_STRING([--disable-matlab], [disable compilation of MEX files for MATLAB]), [], [enable_matlab=yes]) AC_ARG_ENABLE([matlab], AS_HELP_STRING([--disable-matlab], [disable compilation of MEX files for MATLAB]), [], [enable_matlab=yes])
......
doc/manual.xml
View file @ da322f81
...@@ -3610,23 +3610,13 @@ oo_.posterior_hpdsup.measurement_errors_corr.gdp_conso ...@@ -3610,23 +3610,13 @@ oo_.posterior_hpdsup.measurement_errors_corr.gdp_conso
</cmdsynopsis> </cmdsynopsis>
</refsynopsisdiv> </refsynopsisdiv>
<refsect1><title>Options</title> <refsect1><title>Option</title>
<variablelist> <variablelist>
<varlistentry> <varlistentry>
<term><option>parameter_set</option> = <option>prior_mode</option> | <option>prior_mean</option> | <option>posterior_mode</option> | <term><option>parameter_set</option> = <option>prior_mode</option> | <option>prior_mean</option> | <option>posterior_mode</option> |
<option>posterior_mean</option> | <option>posterior_median</option></term> <option>posterior_mean</option> | <option>posterior_median</option></term>
<listitem><para>Specify the parameter set to use for running the smoother. Default value: <option>posterior_mean</option> if Metropolis has been run, else <option>posterior_mode</option>.</para></listitem> <listitem><para>Specify the parameter set to use for running the smoother. Default value: <option>posterior_mean</option> if Metropolis has been run, else <option>posterior_mode</option>.</para></listitem>
</varlistentry> </varlistentry>
<varlistentry>
<term><option>shocks</option> = (
<replaceable>VARIABLE_NAME</replaceable> [<replaceable>VARIABLE_NAME</replaceable> ...]
[ ; <replaceable>VARIABLE_NAME</replaceable> [<replaceable>VARIABLE_NAME</replaceable> ...] ...] )</term>
<listitem><para>...</para></listitem>
</varlistentry>
<varlistentry>
<term><option>labels</option> = ( <replaceable>VARIABLE_NAME</replaceable> [<replaceable>VARIABLE_NAME</replaceable> ...] )</term>
<listitem><para>...</para></listitem>
</varlistentry>
</variablelist> </variablelist>
</refsect1> </refsect1>
...@@ -3884,7 +3874,7 @@ plot_conditional_forecast(periods = 10) e u; ...@@ -3884,7 +3874,7 @@ plot_conditional_forecast(periods = 10) e u;
<refsynopsisdiv> <refsynopsisdiv>
<cmdsynopsis> <cmdsynopsis>
<command>plot_conditional_forecast</command> <command>plot_conditional_forecast</command>
<arg>(<option>periods</option> = <replaceable>INTEGER</replaceable>)</arg><arg choice="plain">;</arg> <arg>(<option>periods</option> = <replaceable>INTEGER</replaceable>)</arg> <arg rep="repeat"><replaceable>VARIABLE_NAME</replaceable></arg><arg choice="plain">;</arg>
</cmdsynopsis> </cmdsynopsis>
</refsynopsisdiv> </refsynopsisdiv>
......
license.txt
View file @ da322f81
...@@ -322,20 +322,3 @@ License: GPL-3+ ...@@ -322,20 +322,3 @@ License: GPL-3+
You should have received a copy of the GNU General Public License You should have received a copy of the GNU General Public License
along with this program. If not, see along with this program. If not, see
<http://www.gnu.org/licenses/>. <http://www.gnu.org/licenses/>.
Files: mex/sources/libslicot/*
Copyright: 2002-2009, NICONET e.V.
License: GPL-2+
This program 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 2 of
the License, or (at your option) any later version.
.
This program 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 this program. If not, see
<http://www.gnu.org/licenses/>.
matlab/DsgeLikelihood.m
View file @ da322f81
...@@ -206,7 +206,7 @@ if (kalman_algo==1)% Multivariate Kalman Filter ...@@ -206,7 +206,7 @@ if (kalman_algo==1)% Multivariate Kalman Filter
end end
if (kalman_algo==2)% Univariate Kalman Filter if (kalman_algo==2)% Univariate Kalman Filter
no_correlation_flag = 1; no_correlation_flag = 1;
if length(H)==1 & H == 0 if isequal(H,0)
H = zeros(nobs,1); H = zeros(nobs,1);
else else
if all(all(abs(H-diag(diag(H)))<1e-14))% ie, the covariance matrix is diagonal... if all(all(abs(H-diag(diag(H)))<1e-14))% ie, the covariance matrix is diagonal...
...@@ -237,7 +237,7 @@ if (kalman_algo==3)% Multivariate Diffuse Kalman Filter ...@@ -237,7 +237,7 @@ if (kalman_algo==3)% Multivariate Diffuse Kalman Filter
end end
if (kalman_algo==4)% Univariate Diffuse Kalman Filter if (kalman_algo==4)% Univariate Diffuse Kalman Filter
no_correlation_flag = 1; no_correlation_flag = 1;
if length(H)==1 & H == 0 if isequal(H,0)
H = zeros(nobs,1); H = zeros(nobs,1);
else else
if all(all(abs(H-diag(diag(H)))<1e-14))% ie, the covariance matrix is diagonal... if all(all(abs(H-diag(diag(H)))<1e-14))% ie, the covariance matrix is diagonal...
......
matlab/DsgeLikelihood_hh.m
View file @ da322f81
...@@ -181,65 +181,7 @@ elseif options_.lik_init == 3 % Diffuse Kalman filter ...@@ -181,65 +181,7 @@ elseif options_.lik_init == 3 % Diffuse Kalman filter
if kalman_algo ~= 4 if kalman_algo ~= 4
kalman_algo = 3; kalman_algo = 3;
end end
[QT,ST] = schur(T); [Z,ST,R1,QT,Pstar,Pinf] = schur_statespace_transformation(mf,T,R,Q,options_.qz_criterium);
e1 = abs(ordeig(ST)) > 2-options_.qz_criterium;
[QT,ST] = ordschur(QT,ST,e1);
k = find(abs(ordeig(ST)) > 2-options_.qz_criterium);
nk = length(k);
nk1 = nk+1;
Pinf = zeros(np,np);
Pinf(1:nk,1:nk) = eye(nk);
Pstar = zeros(np,np);
B = QT'*R*Q*R'*QT;
for i=np:-1:nk+2
if ST(i,i-1) == 0
if i == np
c = zeros(np-nk,1);
else
c = ST(nk1:i,:)*(Pstar(:,i+1:end)*ST(i,i+1:end)')+...
ST(i,i)*ST(nk1:i,i+1:end)*Pstar(i+1:end,i);
end
q = eye(i-nk)-ST(nk1:i,nk1:i)*ST(i,i);
Pstar(nk1:i,i) = q\(B(nk1:i,i)+c);
Pstar(i,nk1:i-1) = Pstar(nk1:i-1,i)';
else
if i == np
c = zeros(np-nk,1);
c1 = zeros(np-nk,1);
else
c = ST(nk1:i,:)*(Pstar(:,i+1:end)*ST(i,i+1:end)')+...
ST(i,i)*ST(nk1:i,i+1:end)*Pstar(i+1:end,i)+...
ST(i,i-1)*ST(nk1:i,i+1:end)*Pstar(i+1:end,i-1);
c1 = ST(nk1:i,:)*(Pstar(:,i+1:end)*ST(i-1,i+1:end)')+...
ST(i-1,i-1)*ST(nk1:i,i+1:end)*Pstar(i+1:end,i-1)+...
ST(i-1,i)*ST(nk1:i,i+1:end)*Pstar(i+1:end,i);
end
q = [eye(i-nk)-ST(nk1:i,nk1:i)*ST(i,i) -ST(nk1:i,nk1:i)*ST(i,i-1);...
-ST(nk1:i,nk1:i)*ST(i-1,i) eye(i-nk)-ST(nk1:i,nk1:i)*ST(i-1,i-1)];
z = q\[B(nk1:i,i)+c;B(nk1:i,i-1)+c1];
Pstar(nk1:i,i) = z(1:(i-nk));
Pstar(nk1:i,i-1) = z(i-nk+1:end);
Pstar(i,nk1:i-1) = Pstar(nk1:i-1,i)';
Pstar(i-1,nk1:i-2) = Pstar(nk1:i-2,i-1)';
i = i - 1;
end
end
if i == nk+2
c = ST(nk+1,:)*(Pstar(:,nk+2:end)*ST(nk1,nk+2:end)')+ST(nk1,nk1)*ST(nk1,nk+2:end)*Pstar(nk+2:end,nk1);
Pstar(nk1,nk1)=(B(nk1,nk1)+c)/(1-ST(nk1,nk1)*ST(nk1,nk1));
end
Z = QT(mf,:);
R1 = QT'*R;
[QQ,RR,EE] = qr(Z*ST(:,1:nk),0);
k = find(abs(diag([RR; zeros(nk-size(Z,1),size(RR,2))])) < 1e-8);
if length(k) > 0
k1 = EE(:,k);
dd =ones(nk,1);
dd(k1) = zeros(length(k1),1);
Pinf(1:nk,1:nk) = diag(dd);
end
end
if kalman_algo == 2
end end
kalman_tol = options_.kalman_tol; kalman_tol = options_.kalman_tol;
riccati_tol = options_.riccati_tol; riccati_tol = options_.riccati_tol;
...@@ -262,7 +204,7 @@ if (kalman_algo==1)% Multivariate Kalman Filter ...@@ -262,7 +204,7 @@ if (kalman_algo==1)% Multivariate Kalman Filter
end end
if (kalman_algo==2)% Univariate Kalman Filter if (kalman_algo==2)% Univariate Kalman Filter
no_correlation_flag = 1; no_correlation_flag = 1;
if length(H)==1 & H == 0 if isequal(H,0)
H = zeros(nobs,1); H = zeros(nobs,1);
else else
if all(all(abs(H-diag(diag(H)))<1e-14))% ie, the covariance matrix is diagonal... if all(all(abs(H-diag(diag(H)))<1e-14))% ie, the covariance matrix is diagonal...
...@@ -293,7 +235,7 @@ if (kalman_algo==3)% Multivariate Diffuse Kalman Filter ...@@ -293,7 +235,7 @@ if (kalman_algo==3)% Multivariate Diffuse Kalman Filter
end end
if (kalman_algo==4)% Univariate Diffuse Kalman Filter if (kalman_algo==4)% Univariate Diffuse Kalman Filter
no_correlation_flag = 1; no_correlation_flag = 1;
if length(H)==1 & H == 0 if isequal(H,0)
H = zeros(nobs,1); H = zeros(nobs,1);
else else
if all(all(abs(H-diag(diag(H)))<1e-14))% ie, the covariance matrix is diagonal... if all(all(abs(H-diag(diag(H)))<1e-14))% ie, the covariance matrix is diagonal...
......
matlab/DsgeSmoother.m
View file @ da322f81
...@@ -133,7 +133,7 @@ data1 = Y-trend; ...@@ -133,7 +133,7 @@ data1 = Y-trend;
% ----------------------------------------------------------------------------- % -----------------------------------------------------------------------------
% 4. Kalman smoother % 4. Kalman smoother
% ----------------------------------------------------------------------------- % -----------------------------------------------------------------------------
if size(H,1) == 1 && H == 0 if isequal(H,0)
H = zeros(nobs,nobs); H = zeros(nobs,nobs);
end end
...@@ -156,7 +156,7 @@ if kalman_algo == 1 || kalman_algo == 3 ...@@ -156,7 +156,7 @@ if kalman_algo == 1 || kalman_algo == 3
[alphahat,epsilonhat,etahat,ahat,P,aK,PK,d,decomp] = missing_DiffuseKalmanSmootherH1_Z(ST, ... [alphahat,epsilonhat,etahat,ahat,P,aK,PK,d,decomp] = missing_DiffuseKalmanSmootherH1_Z(ST, ...
Z,R1,Q,H,Pinf,Pstar, ... Z,R1,Q,H,Pinf,Pstar, ...
data1,nobs,np,smpl,kalman_tol,riccati_tol,data_index); data1,nobs,np,smpl,kalman_tol,riccati_tol,data_index);
if all(alphahat(:)==0) if isequal(alphahat,0)
if kalman_algo == 1 if kalman_algo == 1
kalman_algo = 2; kalman_algo = 2;
elseif kalman_algo == 3 elseif kalman_algo == 3
...@@ -214,147 +214,7 @@ if estim_params_.ncn && (kalman_algo == 2 || kalman_algo == 4) ...@@ -214,147 +214,7 @@ if estim_params_.ncn && (kalman_algo == 2 || kalman_algo == 4)
if ~isempty(decomp) if ~isempty(decomp)
decomp = decomp(:,k,:,:); decomp = decomp(:,k,:,:);
end end
if ~isempte(P) if ~isempty(P)
P = P(k,k,:); P = P(k,k,:);
end end
end end
% $$$ if any(any(H ~= 0)) % should be replaced by a flag
% $$$ if kalman_algo == 1
% $$$ [alphahat,epsilonhat,etahat,ahat,P,aK,PK,decomp] = ...
% $$$ kalman_smoother(T,R,Q,H,Pstar,data1,start,mf,kalman_tol,riccati_tol);
% $$$ if all(alphahat(:)==0)
% $$$ kalman_algo = 2;
% $$$ if ~estim_params_.ncn
% $$$ [alphahat,epsilonhat,etahat,ahat,aK] = ...
% $$$ DiffuseKalmanSmootherH3(T,R,Q,H,Pinf,Pstar,Y,trend,nobs,np,smpl,mf);
% $$$ else
% $$$ [alphahat,epsilonhat,etahat,ahat,aK] = ...
% $$$ DiffuseKalmanSmootherH3corr(T,R,Q,H,Pinf,Pstar,Y,trend, ...
% $$$ nobs,np,smpl,mf);
% $$$ end
% $$$ end
% $$$ elseif kalman_algo == 2
% $$$ if ~estim_params_.ncn
% $$$ [alphahat,epsilonhat,etahat,ahat,aK] = ...
% $$$ DiffuseKalmanSmootherH3(T,R,Q,H,Pinf,Pstar,Y,trend,nobs,np,smpl,mf);
% $$$ else
% $$$ [alphahat,epsilonhat,etahat,ahat,aK] = ...
% $$$ DiffuseKalmanSmootherH3corr(T,R,Q,H,Pinf,Pstar,Y,trend, ...
% $$$ nobs,np,smpl,mf);
% $$$ end
% $$$ elseif kalman_algo == 3 | kalman_algo == 4
% $$$ data1 = Y - trend;
% $$$ if kalman_algo == 3
% $$$ [alphahat,epsilonhat,etahat,ahat,P,aK,PK,d,decomp] = ...
% $$$ DiffuseKalmanSmootherH1_Z(ST,Z,R1,Q,H,Pinf,Pstar,data1,nobs,np,smpl);
% $$$ if all(alphahat(:)==0)
% $$$ kalman_algo = 4;
% $$$ if ~estim_params_.ncn
% $$$ [alphahat,epsilonhat,etahat,ahat,P,aK,PK,d,decomp] = ...
% $$$ DiffuseKalmanSmootherH3_Z(ST,Z,R1,Q,H,Pinf,Pstar,data1,nobs,np,smpl);
% $$$ else
% $$$ [alphahat,epsilonhat,etahat,ahat,P,aK,PK,d,decomp] = ...
% $$$ DiffuseKalmanSmootherH3corr_Z(ST,Z,R1,Q,H,Pinf,Pstar,data1, ...
% $$$ nobs,np,smpl);
% $$$ end
% $$$ end
% $$$ else
% $$$ if ~estim_params_.ncn
% $$$ [alphahat,epsilonhat,etahat,ahat,P,aK,PK,d,decomp] = ...
% $$$ DiffuseKalmanSmootherH3_Z(ST,Z,R1,Q,H,Pinf,Pstar,data1, ...
% $$$ nobs,np,smpl);
% $$$ else
% $$$ [alphahat,epsilonhat,etahat,ahat,P,aK,PK,d,decomp] = ...
% $$$ DiffuseKalmanSmootherH3corr_Z(ST,Z,R1,Q,H,Pinf,Pstar,data1, ...
% $$$ nobs,np,smpl);
% $$$ end
% $$$ end
% $$$ alphahat = QT*alphahat;
% $$$ ahat = QT*ahat;
% $$$ nk = options_.nk;
% $$$ for jnk=1:nk
% $$$ aK(jnk,:,:) = QT*squeeze(aK(jnk,:,:));
% $$$ for i=1:size(PK,4)
% $$$ PK(jnk,:,:,i) = QT*squeeze(PK(jnk,:,:,i))*QT';
% $$$ end
% $$$ for i=1:size(decomp,4)
% $$$ decomp(jnk,:,:,i) = QT*squeeze(decomp(jnk,:,:,i));
% $$$ end
% $$$ end
% $$$ for i=1:size(P,4)
% $$$ P(:,:,i) = QT*squeeze(P(:,:,i))*QT';
% $$$ end
% $$$ end
% $$$ else
% $$$ H = 0;
% $$$ if kalman_algo == 1
% $$$ if missing_value
% $$$ [alphahat,etahat,ahat,aK] = missing_DiffuseKalmanSmoother1(T,R,Q, ...
% $$$ Pinf,Pstar,Y,trend,nobs,np,smpl,mf,data_index);
% $$$ else
% $$$ [alphahat,epsilonhat,etahat,ahat,P,aK,PK,decomp] = ...
% $$$ kalman_smoother(T,R,Q,H,Pstar,data1,start,mf,kalman_tol,riccati_tol);
% $$$ end
% $$$ if all(alphahat(:)==0)
% $$$ kalman_algo = 2;
% $$$ end
% $$$ end
% $$$ if kalman_algo == 2
% $$$ if missing_value
% $$$ [alphahat,etahat,ahat,P,aK,PK,d,decomp] = missing_DiffuseKalmanSmoother3(T,R,Q, ...
% $$$ Pinf,Pstar,Y,trend,nobs,np,smpl,mf,data_index);
% $$$ else
% $$$ [alphahat,etahat,ahat,P,aK,PK,d,decomp] = DiffuseKalmanSmoother3(T,R,Q, ...
% $$$ Pinf,Pstar,Y,trend,nobs,np,smpl,mf);
% $$$ end
% $$$ end
% $$$ if kalman_algo == 3
% $$$ data1 = Y - trend;
% $$$ if missing_value
% $$$ [alphahat,etahat,ahat,P,aK,PK,d,decomp] = missing_DiffuseKalmanSmoother1_Z(ST, ...
% $$$ Z,R1,Q,Pinf,Pstar,data1,nobs,np,smpl,data_index);
% $$$ else
% $$$ [alphahat,etahat,ahat,P,aK,PK,d,decomp] = DiffuseKalmanSmoother1_Z(ST, ...
% $$$ Z,R1,Q,Pinf,Pstar, ...
% $$$ data1,nobs,np,smpl);
% $$$ end
% $$$ if all(alphahat(:)==0)
% $$$ kalman_algo = 4;
% $$$ end
% $$$ end
% $$$ if kalman_algo == 4
% $$$ data1 = Y - trend;
% $$$ if missing_value
% $$$ [alphahat,etahat,ahat,P,aK,PK,d,decomp] = missing_DiffuseKalmanSmoother3_Z(ST, ...
% $$$ Z,R1,Q,Pinf,Pstar,data1,nobs,np,smpl,data_index);
% $$$ else
% $$$ [alphahat,etahat,ahat,P,aK,PK,d,decomp] = DiffuseKalmanSmoother3_Z(ST, ...
% $$$ Z,R1,Q,Pinf,Pstar, ...
% $$$ data1,nobs,np,smpl);
% $$$ end
% $$$ end
% $$$ if kalman_algo == 3 | kalman_algo == 4
% $$$ alphahat = QT*alphahat;
% $$$ ahat = QT*ahat;
% $$$ nk = options_.nk;
% $$$ % $$$ if M_.exo_nbr<2 % Fix the crash of Dynare when the estimated model has only one structural shock (problem with
% $$$ % $$$ % the squeeze function, that does not affect 2D arrays).
% $$$ % $$$ size_decomp = 0;
% $$$ % $$$ else
% $$$ % $$$ size_decomp = size(decomp,4);
% $$$ % $$$ end
% $$$ for jnk=1:nk
% $$$ aK(jnk,:,:) = QT*squeeze(aK(jnk,:,:));
% $$$ for i=1:size(PK,4)
% $$$ PK(jnk,:,:,i) = QT*dynare_squeeze(PK(jnk,:,:,i))*QT';
% $$$ end
% $$$ for i=1:size(decomp,4)
% $$$ decomp(jnk,:,:,i) = QT*dynare_squeeze(decomp(jnk,:,:,i));
% $$$ end
% $$$ end
% $$$ for i=1:size(P,4)
% $$$ P(:,:,i) = QT*dynare_squeeze(P(:,:,i))*QT';
% $$$ end
% $$$ end
% $$$ end
matlab/dynare_config.m
View file @ da322f81
...@@ -43,7 +43,6 @@ addpath([dynareroot '/kalman/likelihood']) ...@@ -43,7 +43,6 @@ addpath([dynareroot '/kalman/likelihood'])
addpath([dynareroot '/kalman/smoother']) addpath([dynareroot '/kalman/smoother'])
addpath([dynareroot '/AIM/']) addpath([dynareroot '/AIM/'])
addpath([dynareroot '/partial_information/']) addpath([dynareroot '/partial_information/'])
addpath([dynareroot '/ms-sbvar/'])
addpath([dynareroot '/parallel/']) addpath([dynareroot '/parallel/'])
% For functions that exist only under some Octave versions % For functions that exist only under some Octave versions
......
matlab/dynare_estimation_1.m
View file @ da322f81
...@@ -220,7 +220,7 @@ else ...@@ -220,7 +220,7 @@ else
k3 = (1:M_.endo_nbr)'; k3 = (1:M_.endo_nbr)';
end end
bayestopt_.smoother_var_list = union(k2,k3); bayestopt_.smoother_var_list = union(k2,k3);
[junk,bayestopt_.smoother_saved_var_list] = intersect(k3,bayestopt_.smoother_var_list); [junk,bayestopt_.smoother_saved_var_list] = intersect(k3,bayestopt_.smoother_var_list(:));
[junk,ic] = intersect(bayestopt_.smoother_var_list,nstatic+(1:npred)'); [junk,ic] = intersect(bayestopt_.smoother_var_list,nstatic+(1:npred)');
bayestopt_.smoother_restrict_columns = ic; bayestopt_.smoother_restrict_columns = ic;
[junk,bayestopt_.smoother_mf] = ismember(var_obs_index, ... [junk,bayestopt_.smoother_mf] = ismember(var_obs_index, ...
...@@ -306,7 +306,7 @@ options_ = set_default_option(options_,'mh_nblck',2); ...@@ -306,7 +306,7 @@ options_ = set_default_option(options_,'mh_nblck',2);
options_ = set_default_option(options_,'nodiagnostic',0); options_ = set_default_option(options_,'nodiagnostic',0);
% load mode file is necessary % load mode file is necessary
if length(options_.mode_file) > 0 && ~options_.mh_posterior_mode_estimation if ~isempty(options_.mode_file) && ~options_.mh_posterior_mode_estimation
load(options_.mode_file); load(options_.mode_file);
end end
...@@ -366,7 +366,7 @@ missing_value = ~(number_of_observations == gend*n_varobs); ...@@ -366,7 +366,7 @@ missing_value = ~(number_of_observations == gend*n_varobs);
initial_estimation_checks(xparam1,gend,data,data_index,number_of_observations,no_more_missing_observations); initial_estimation_checks(xparam1,gend,data,data_index,number_of_observations,no_more_missing_observations);
if isnumeric(options_.mode_compute) && options_.mode_compute == 0 && length(options_.mode_file) == 0 && options_.mh_posterior_mode_estimation==0 if isequal(options_.mode_compute,0) && isempty(options_.mode_file) && options_.mh_posterior_mode_estimation==0
if options_.smoother == 1 if options_.smoother == 1
[atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,T,R,P,PK,decomp] = DsgeSmoother(xparam1,gend,data,data_index,missing_value); [atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,T,R,P,PK,decomp] = DsgeSmoother(xparam1,gend,data,data_index,missing_value);
oo_.Smoother.SteadyState = ys; oo_.Smoother.SteadyState = ys;
...@@ -401,14 +401,14 @@ if isnumeric(options_.mode_compute) && options_.mode_compute == 0 && length(opti ...@@ -401,14 +401,14 @@ if isnumeric(options_.mode_compute) && options_.mode_compute == 0 && length(opti
return; return;
end end
if options_.mode_compute==6 if isequal(options_.mode_compute,6)
% Erase previously computed optimal mh scale parameter. % Erase previously computed optimal mh scale parameter.
delete([M_.fname '_optimal_mh_scale_parameter.mat']) delete([M_.fname '_optimal_mh_scale_parameter.mat'])
end end
%% Estimation of the posterior mode or likelihood mode %% Estimation of the posterior mode or likelihood mode
if any(options_.mode_compute ~= 0) && ~options_.mh_posterior_mode_estimation if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation
if ~options_.dsge_var if ~options_.dsge_var
fh=str2func('DsgeLikelihood'); fh=str2func('DsgeLikelihood');
else else
...@@ -470,7 +470,7 @@ if any(options_.mode_compute ~= 0) && ~options_.mh_posterior_mode_estimation ...@@ -470,7 +470,7 @@ if any(options_.mode_compute ~= 0) && ~options_.mh_posterior_mode_estimation
if isfield(options_,'ftol') if isfield(options_,'ftol')
crit = options_.ftol; crit = options_.ftol;
else else
crit = 1.e-7; crit = 1.e-5;
end end
if isfield(options_,'nit') if isfield(options_,'nit')
nit = options_.nit; nit = options_.nit;
...@@ -645,8 +645,7 @@ if any(options_.mode_compute ~= 0) && ~options_.mh_posterior_mode_estimation ...@@ -645,8 +645,7 @@ if any(options_.mode_compute ~= 0) && ~options_.mh_posterior_mode_estimation
' option is unknown!']) ' option is unknown!'])
end end
end end
% if options_.mode_compute ~= 5 if ~isequal(options_.mode_compute,6)
if options_.mode_compute ~= 6
if options_.cova_compute == 1 if options_.cova_compute == 1
if ~options_.dsge_var if ~options_.dsge_var
hh = reshape(hessian('DsgeLikelihood',xparam1, ... hh = reshape(hessian('DsgeLikelihood',xparam1, ...
...@@ -665,7 +664,6 @@ if any(options_.mode_compute ~= 0) && ~options_.mh_posterior_mode_estimation ...@@ -665,7 +664,6 @@ if any(options_.mode_compute ~= 0) && ~options_.mh_posterior_mode_estimation
else else
save([M_.fname '_mode.mat'],'xparam1','parameter_names'); save([M_.fname '_mode.mat'],'xparam1','parameter_names');
end end
% end
end end
if options_.cova_compute == 0 if options_.cova_compute == 0
...@@ -911,7 +909,7 @@ if any(bayestopt_.pshape > 0) & options_.TeX %% Bayesian estimation (posterior m ...@@ -911,7 +909,7 @@ if any(bayestopt_.pshape > 0) & options_.TeX %% Bayesian estimation (posterior m
fprintf(fidTeX,['%% ' datestr(now,0)]); fprintf(fidTeX,['%% ' datestr(now,0)]);
fprintf(fidTeX,' \n'); fprintf(fidTeX,' \n');
fprintf(fidTeX,' \n'); fprintf(fidTeX,' \n');
fprintf(fidTeX,'{\\tiny \n') fprintf(fidTeX,'{\\tiny \n');
fprintf(fidTeX,'\\begin{table}\n'); fprintf(fidTeX,'\\begin{table}\n');
fprintf(fidTeX,'\\centering\n'); fprintf(fidTeX,'\\centering\n');
fprintf(fidTeX,'\\begin{tabular}{l|lcccc} \n'); fprintf(fidTeX,'\\begin{tabular}{l|lcccc} \n');
...@@ -934,7 +932,7 @@ if any(bayestopt_.pshape > 0) & options_.TeX %% Bayesian estimation (posterior m ...@@ -934,7 +932,7 @@ if any(bayestopt_.pshape > 0) & options_.TeX %% Bayesian estimation (posterior m
fprintf(fidTeX,'\\caption{Results from posterior parameters (parameters)}\n '); fprintf(fidTeX,'\\caption{Results from posterior parameters (parameters)}\n ');
fprintf(fidTeX,'\\label{Table:Posterior:1}\n'); fprintf(fidTeX,'\\label{Table:Posterior:1}\n');
fprintf(fidTeX,'\\end{table}\n'); fprintf(fidTeX,'\\end{table}\n');
fprintf(fidTeX,'} \n') fprintf(fidTeX,'} \n');
fprintf(fidTeX,'%% End of TeX file.\n'); fprintf(fidTeX,'%% End of TeX file.\n');
fclose(fidTeX); fclose(fidTeX);
end end
...@@ -951,7 +949,7 @@ if any(bayestopt_.pshape > 0) & options_.TeX %% Bayesian estimation (posterior m ...@@ -951,7 +949,7 @@ if any(bayestopt_.pshape > 0) & options_.TeX %% Bayesian estimation (posterior m
fprintf(fidTeX,'\\centering\n'); fprintf(fidTeX,'\\centering\n');
fprintf(fidTeX,'\\begin{tabular}{l|lcccc} \n'); fprintf(fidTeX,'\\begin{tabular}{l|lcccc} \n');
fprintf(fidTeX,'\\hline\\hline \\\\ \n'); fprintf(fidTeX,'\\hline\\hline \\\\ \n');
fprintf(fidTeX,' & Prior distribution & Prior mean & Prior s.d. & Posterior mode & s.d. \\\\ \n') fprintf(fidTeX,' & Prior distribution & Prior mean & Prior s.d. & Posterior mode & s.d. \\\\ \n');
fprintf(fidTeX,'\\hline \\\\ \n'); fprintf(fidTeX,'\\hline \\\\ \n');
ip = 1; ip = 1;
for i=1:nvx for i=1:nvx
...@@ -970,7 +968,7 @@ if any(bayestopt_.pshape > 0) & options_.TeX %% Bayesian estimation (posterior m ...@@ -970,7 +968,7 @@ if any(bayestopt_.pshape > 0) & options_.TeX %% Bayesian estimation (posterior m
fprintf(fidTeX,'\\caption{Results from posterior parameters (standard deviation of structural shocks)}\n '); fprintf(fidTeX,'\\caption{Results from posterior parameters (standard deviation of structural shocks)}\n ');
fprintf(fidTeX,'\\label{Table:Posterior:2}\n'); fprintf(fidTeX,'\\label{Table:Posterior:2}\n');
fprintf(fidTeX,'\\end{table}\n'); fprintf(fidTeX,'\\end{table}\n');
fprintf(fidTeX,'} \n') fprintf(fidTeX,'} \n');
fprintf(fidTeX,'%% End of TeX file.\n'); fprintf(fidTeX,'%% End of TeX file.\n');
fclose(fidTeX); fclose(fidTeX);
end end
......
matlab/ms-sbvar/clean_swz_files.m deleted 100644 → 0
View file @ b00b2fad
function clean_swz_files(mod_name)
% function clean_swz_files()
% removes SWZ intermediary files
%
% INPUTS
% none
%
% OUTPUTS
% none
%
% SPECIAL REQUIREMENTS
% none
% Copyright (C) 2007-2008, 2010 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/>.
delete_if_exist(['./draws_' mod_name '.dat'])
delete_if_exist(['est_aux_' mod_name '.dat'])
delete_if_exist(['est_csminwel_' mod_name '.dat'])
delete_if_exist(['est_final_' mod_name '.dat'])
delete_if_exist(['est_flat_header_' mod_name '.dat'])
delete_if_exist(['est_flat_' mod_name '.dat'])
delete_if_exist(['est_intermediate_' mod_name '.dat'])
delete_if_exist(['header_' mod_name '.dat'])
delete_if_exist(['init_' mod_name '.dat'])
delete_if_exist('markov_file.dat')
delete_if_exist(['matlab_' mod_name '.prn'])
delete_if_exist(['mhm_draws_' mod_name '.dat'])
delete_if_exist(['mhm_input_' mod_name '.dat'])
delete_if_exist(['mhm_intermediate_draws_' mod_name '.dat'])
delete_if_exist(['mhm_intermediate_' mod_name '.dat'])
delete_if_exist(['mhm_regime_counts_' mod_name '.dat'])
delete_if_exist(['probabilities_' mod_name '.dat'])
delete_if_exist(['truncatedpower_md_posterior_' mod_name '.dat'])
delete_if_exist(['truncatedpower_md_proposal_' mod_name '.dat'])
delete_if_exist(['truncatedpower_md_' mod_name '.dat'])
function delete_if_exist(fname)
if exist(fname) == 2
delete(fname)
end
\ No newline at end of file
matlab/ms-sbvar/cstz/bfgsi.m deleted 100644 → 0
View file @ b00b2fad
function H = bfgsi(H0,dg,dx)
% H = bfgsi(H0,dg,dx)
% dg is previous change in gradient; dx is previous change in x;
% 6/8/93 version that updates inverse hessian instead of hessian
% itself.
% Copyright by Christopher Sims 1996. This material may be freely
% reproduced and modified.
dispIndx = 0; % 1: turn on all the diplays on the screen; 0: turn off (Added by T. Zha)
if size(dg,2)>1
dg=dg';
end
if size(dx,2)>1
dx=dx';
end
Hdg = H0*dg;
dgdx = dg'*dx;
if (abs(dgdx) >1e-12)
H = H0 + (1+(dg'*Hdg)/dgdx)*(dx*dx')/dgdx - (dx*Hdg'+Hdg*dx')/dgdx;
else
if dispIndx
disp('bfgs update failed.')
disp(['|dg| = ' num2str(sqrt(dg'*dg)) '|dx| = ' num2str(sqrt(dx'*dx))]);
disp(['dg''*dx = ' num2str(dgdx)])
disp(['|H*dg| = ' num2str(Hdg'*Hdg)])
end
H=H0;
end
save H.dat H
matlab/ms-sbvar/cstz/csminit.m deleted 100644 → 0
View file @ b00b2fad
function [fhat,xhat,fcount,retcode] = csminit(fcn,x0,f0,g0,badg,H0,varargin)
% [fhat,xhat,fcount,retcode] = csminit(fcn,x0,f0,g0,badg,H0,...
% P1,P2,P3,P4,P5,P6,P7,P8)
% retcodes: 0, normal step. 5, largest step still improves too fast.
% 4,2 back and forth adjustment of stepsize didn't finish. 3, smallest
% stepsize still improves too slow. 6, no improvement found. 1, zero
% gradient.
%---------------------
% Modified 7/22/96 to omit variable-length P list, for efficiency and compilation.
% Places where the number of P's need to be altered or the code could be returned to
% its old form are marked with ARGLIST comments.
%
% Fixed 7/17/93 to use inverse-hessian instead of hessian itself in bfgs
% update.
%
% Fixed 7/19/93 to flip eigenvalues of H to get better performance when
% it's not psd.
%
% Fixed 02/19/05 to correct for low angle problems.
%
%tailstr = ')';
%for i=nargin-6:-1:1
% tailstr=[ ',P' num2str(i) tailstr];
%end
dispIndx = 0; % 1: turn on all the diplays on the screen; 0: turn off (Added by T. Zha)
%ANGLE = .03; % when output of this variable becomes negative, we have wrong analytical graident
ANGLE = .005; % works for identified VARs and OLS
%THETA = .03;
THETA = .3; %(0<THETA<.5) THETA near .5 makes long line searches, possibly fewer iterations.
FCHANGE = 1000;
MINLAMB = 1e-9;
% fixed 7/15/94
% MINDX = .0001;
% MINDX = 1e-6;
MINDFAC = .01;
fcount=0;
lambda=1;
xhat=x0;
f=f0;
fhat=f0;
g = g0;
gnorm = norm(g);
%
if (gnorm < 1.e-12) & ~badg % put ~badg 8/4/94
retcode =1;
dxnorm=0;
% gradient convergence
else
% with badg true, we don't try to match rate of improvement to directional
% derivative. We're satisfied just to get some improvement in f.
%
%if(badg)
% dx = -g*FCHANGE/(gnorm*gnorm);
% dxnorm = norm(dx);
% if dxnorm > 1e12
% disp('Bad, small gradient problem.')
% dx = dx*FCHANGE/dxnorm;
% end
%else
% Gauss-Newton step;
%---------- Start of 7/19/93 mod ---------------
%[v d] = eig(H0);
%toc
%d=max(1e-10,abs(diag(d)));
%d=abs(diag(d));
%dx = -(v.*(ones(size(v,1),1)*d'))*(v'*g);
% toc
dx = -H0*g;
% toc
dxnorm = norm(dx);
if dxnorm > 1e12
if dispIndx, disp('Near-singular H problem.'), end
dx = dx*FCHANGE/dxnorm;
end
dfhat = dx'*g0;
%end
%
%
if ~badg
% test for alignment of dx with gradient and fix if necessary
a = -dfhat/(gnorm*dxnorm);
if a<ANGLE
dx = dx - (ANGLE*dxnorm/gnorm+dfhat/(gnorm*gnorm))*g;
% suggested alternate code: ---------------------
dx = dx*dxnorm/norm(dx); % Added 02/19/05 by CAS. This keeps scale invariant to the angle correction
% ------------------------------------------------
dfhat = dx'*g;
% dxnorm = norm(dx); % Removed 02/19/05 by CAS. This line unnecessary with modification that keeps scale invariant
if dispIndx, disp(sprintf('Correct for low angle: %g',a)), end
end
end
if dispIndx, disp(sprintf('Predicted improvement: %18.9f',-dfhat/2)), end
%
% Have OK dx, now adjust length of step (lambda) until min and
% max improvement rate criteria are met.
done=0;
factor=3;
shrink=1;
lambdaMin=0;
lambdaMax=inf;
lambdaPeak=0;
fPeak=f0;
lambdahat=0;
while ~done
if size(x0,2)>1
dxtest=x0+dx'*lambda;
else
dxtest=x0+dx*lambda;
end
% home
f = feval(fcn,dxtest,varargin{:});
%ARGLIST
%f = feval(fcn,dxtest,P1,P2,P3,P4,P5,P6,P7,P8,P9,P10,P11,P12,P13);
% f = feval(fcn,x0+dx*lambda,P1,P2,P3,P4,P5,P6,P7,P8);
if dispIndx, disp(sprintf('lambda = %10.5g; f = %20.7f',lambda,f )), end
%debug
%disp(sprintf('Improvement too great? f0-f: %g, criterion: %g',f0-f,-(1-THETA)*dfhat*lambda))
if f<fhat
fhat=f;
xhat=dxtest;
lambdahat = lambda;
end
fcount=fcount+1;
shrinkSignal = (~badg & (f0-f < max([-THETA*dfhat*lambda 0]))) | (badg & (f0-f) < 0) ;
growSignal = ~badg & ( (lambda > 0) & (f0-f > -(1-THETA)*dfhat*lambda) );
if shrinkSignal & ( (lambda>lambdaPeak) | (lambda<0) )
if (lambda>0) & ((~shrink) | (lambda/factor <= lambdaPeak))
shrink=1;
factor=factor^.6;
while lambda/factor <= lambdaPeak
factor=factor^.6;
end
%if (abs(lambda)*(factor-1)*dxnorm < MINDX) | (abs(lambda)*(factor-1) < MINLAMB)
if abs(factor-1)<MINDFAC
if abs(lambda)<4
retcode=2;
else
retcode=7;
end
done=1;
end
end
if (lambda<lambdaMax) & (lambda>lambdaPeak)
lambdaMax=lambda;
end
lambda=lambda/factor;
if abs(lambda) < MINLAMB
if (lambda > 0) & (f0 <= fhat)
% try going against gradient, which may be inaccurate
if dispIndx, lambda = -lambda*factor^6, end
else
if lambda < 0
retcode = 6;
else
retcode = 3;
end
done = 1;
end
end
elseif (growSignal & lambda>0) | (shrinkSignal & ((lambda <= lambdaPeak) & (lambda>0)))
if shrink
shrink=0;
factor = factor^.6;
%if ( abs(lambda)*(factor-1)*dxnorm< MINDX ) | ( abs(lambda)*(factor-1)< MINLAMB)
if abs(factor-1)<MINDFAC
if abs(lambda)<4
retcode=4;
else
retcode=7;
end
done=1;
end
end
if ( f<fPeak ) & (lambda>0)
fPeak=f;
lambdaPeak=lambda;
if lambdaMax<=lambdaPeak
lambdaMax=lambdaPeak*factor*factor;
end
end
lambda=lambda*factor;
if abs(lambda) > 1e20;
retcode = 5;
done =1;
end
else
done=1;
if factor < 1.2
retcode=7;
else
retcode=0;
end
end
end
end
if dispIndx, disp(sprintf('Norm of dx %10.5g', dxnorm)), end
matlab/ms-sbvar/cstz/csminwel.m deleted 100644 → 0
View file @ b00b2fad
function [fh,xh,gh,H,itct,fcount,retcodeh] = csminwel(fcn,x0,H0,grad,crit,nit,varargin)
%[fhat,xhat,ghat,Hhat,itct,fcount,retcodehat] = csminwel(fcn,x0,H0,grad,crit,nit,varargin)
% fcn: string naming the objective function to be minimized
% x0: initial value of the parameter vector
% H0: initial value for the inverse Hessian. Must be positive definite.
% grad: Either a string naming a function that calculates the gradient, or the null matrix.
% If it's null, the program calculates a numerical gradient. In this case fcn must
% be written so that it can take a matrix argument and produce a row vector of values.
% crit: Convergence criterion. Iteration will cease when it proves impossible to improve the
% function value by more than crit.
% nit: Maximum number of iterations.
% varargin: A list of optional length of additional parameters that get handed off to fcn each
% time it is called.
% Note that if the program ends abnormally, it is possible to retrieve the current x,
% f, and H from the files g1.mat and H.mat that are written at each iteration and at each
% hessian update, respectively. (When the routine hits certain kinds of difficulty, it
% write g2.mat and g3.mat as well. If all were written at about the same time, any of them
% may be a decent starting point. One can also start from the one with best function value.)
% NOTE: The display on screen can be turned off by seeting dispIndx=0 in this
% function. This option is used for the loop operation. T. Zha, 2 May 2000
% NOTE: You may want to change stps to 1.0e-02 or 1.0e-03 to get a better convergence. August, 2006
Verbose = 0; % 1: turn on all the diplays on the screen; 0: turn off (Added by T. Zha)
dispIndx = 0; % 1: turn on all the diplays on the screen; 0: turn off (Added by T. Zha)
[nx,no]=size(x0);
nx=max(nx,no);
NumGrad= ( ~ischar(grad) | length(grad)==0);
done=0;
itct=0;
fcount=0;
snit=100;
%tailstr = ')';
%stailstr = [];
% Lines below make the number of Pi's optional. This is inefficient, though, and precludes
% use of the matlab compiler. Without them, we use feval and the number of Pi's must be
% changed with the editor for each application. Places where this is required are marked
% with ARGLIST comments
%for i=nargin-6:-1:1
% tailstr=[ ',P' num2str(i) tailstr];
% stailstr=[' P' num2str(i) stailstr];
%end
f0 = eval([fcn '(x0,varargin{:})']);
%ARGLIST
%f0 = feval(fcn,x0,P1,P2,P3,P4,P5,P6,P7,P8,P9,P10,P11,P12,P13);
% disp('first fcn in csminwel.m ----------------') % Jinill on 9/5/95
if f0 > 1e50, disp('Bad initial parameter.'), return, end
if NumGrad
if length(grad)==0
[g badg] = numgradcd(fcn,x0, varargin{:});
%ARGLIST
%[g badg] = numgradcd(fcn,x0,P1,P2,P3,P4,P5,P6,P7,P8,P9,P10,P11,P12,P13);
else
badg=any(find(grad==0));
g=grad;
end
%numgradcd(fcn,x0,P1,P2,P3,P4);
else
[g badg] = eval([grad '(x0,varargin{:})']);
%ARGLIST
%[g badg] = feval(grad,x0,P1,P2,P3,P4,P5,P6,P7,P8,P9,P10,P11,P12,P13);
end
retcode3=101;
x=x0;
f=f0;
H=H0;
cliff=0;
while ~done
g1=[]; g2=[]; g3=[];
%addition fj. 7/6/94 for control
if dispIndx
disp('-----------------')
disp('-----------------')
%disp('f and x at the beginning of new iteration')
disp(sprintf('f at the beginning of new iteration, %20.10f',f))
%-----------Comment out this line if the x vector is long----------------
disp([sprintf('x = ') sprintf('%15.8g%15.8g%15.8g%15.8g%15.8g\n',x)]);
end
%-------------------------
itct=itct+1;
[f1 x1 fc retcode1] = csminit(fcn,x,f,g,badg,H,varargin{:});
%ARGLIST
%[f1 x1 fc retcode1] = csminit(fcn,x,f,g,badg,H,P1,P2,P3,P4,P5,P6,P7,...
% P8,P9,P10,P11,P12,P13);
% itct=itct+1;
fcount = fcount+fc;
% erased on 8/4/94
% if (retcode == 1) | (abs(f1-f) < crit)
% done=1;
% end
% if itct > nit
% done = 1;
% retcode = -retcode;
% end
if retcode1 ~= 1
if retcode1==2 | retcode1==4
wall1=1; badg1=1;
else
if NumGrad
[g1 badg1] = numgradcd(fcn, x1,varargin{:});
%ARGLIST
%[g1 badg1] = numgradcd(fcn, x1,P1,P2,P3,P4,P5,P6,P7,P8,P9,...
% P10,P11,P12,P13);
else
[g1 badg1] = eval([grad '(x1,varargin{:})']);
%ARGLIST
%[g1 badg1] = feval(grad, x1,P1,P2,P3,P4,P5,P6,P7,P8,P9,...
% P10,P11,P12,P13);
end
wall1=badg1;
% g1
save g1 g1 x1 f1 varargin;
%ARGLIST
%save g1 g1 x1 f1 P1 P2 P3 P4 P5 P6 P7 P8 P9 P10 P11 P12 P13;
end
if wall1 % & (~done) by Jinill
% Bad gradient or back and forth on step length. Possibly at
% cliff edge. Try perturbing search direction.
%
%fcliff=fh;xcliff=xh;
if dispIndx
disp(' ')
disp('************************* Random search. *****************************************')
disp('************************* Random search. *****************************************')
disp(' ')
pause(1.0)
end
Hcliff=H+diag(diag(H).*rand(nx,1));
if dispIndx, disp('Cliff. Perturbing search direction.'), end
[f2 x2 fc retcode2] = csminit(fcn,x,f,g,badg,Hcliff,varargin{:});
%ARGLIST
%[f2 x2 fc retcode2] = csminit(fcn,x,f,g,badg,Hcliff,P1,P2,P3,P4,...
% P5,P6,P7,P8,P9,P10,P11,P12,P13);
fcount = fcount+fc; % put by Jinill
if f2 < f
if retcode2==2 | retcode2==4
wall2=1; badg2=1;
else
if NumGrad
[g2 badg2] = numgradcd(fcn, x2,varargin{:});
%ARGLIST
%[g2 badg2] = numgradcd(fcn, x2,P1,P2,P3,P4,P5,P6,P7,P8,...
% P9,P10,P11,P12,P13);
else
[g2 badg2] = eval([grad '(x2,varargin{:})']);
%ARGLIST
%[g2 badg2] = feval(grad,x2,P1,P2,P3,P4,P5,P6,P7,P8,...
% P9,P10,P11,P12,P13);
end
wall2=badg2;
% g2
if dispIndx, badg2, end
save g2 g2 x2 f2 varargin
%ARGLIST
%save g2 g2 x2 f2 P1 P2 P3 P4 P5 P6 P7 P8 P9 P10 P11 P12 P13;
end
if wall2
if dispIndx, disp('Cliff again. Try traversing'), end
if norm(x2-x1) < 1e-13
f3=f; x3=x; badg3=1;retcode3=101;
else
gcliff=((f2-f1)/((norm(x2-x1))^2))*(x2-x1);
[f3 x3 fc retcode3] = csminit(fcn,x,f,gcliff,0,eye(nx),varargin{:});
%ARGLIST
%[f3 x3 fc retcode3] = csminit(fcn,x,f,gcliff,0,eye(nx),P1,P2,P3,...
% P4,P5,P6,P7,P8,...
% P9,P10,P11,P12,P13);
fcount = fcount+fc; % put by Jinill
if retcode3==2 | retcode3==4
wall3=1; badg3=1;
else
if NumGrad
[g3 badg3] = numgradcd(fcn, x3,varargin{:});
%ARGLIST
%[g3 badg3] = numgradcd(fcn, x3,P1,P2,P3,P4,P5,P6,P7,P8,...
% P9,P10,P11,P12,P13);
else
[g3 badg3] = eval([grad '(x3,varargin{:})']);
%ARGLIST
%[g3 badg3] = feval(grad,x3,P1,P2,P3,P4,P5,P6,P7,P8,...
% P9,P10,P11,P12,P13);
end
wall3=badg3;
% g3
if dispIndx, badg3, end
save g3 g3 x3 f3 varargin;
%ARGLIST
%save g3 g3 x3 f3 P1 P2 P3 P4 P5 P6 P7 P8 P9 P10 P11 P12 P13;
end
end
else
f3=f; x3=x; badg3=1; retcode3=101;
end
else
f3=f; x3=x; badg3=1;retcode3=101;
end
else
% normal iteration, no walls, or else we're finished here.
f2=f; f3=f; badg2=1; badg3=1; retcode2=101; retcode3=101;
end
else
f1=f; f2=f; f3=f; retcode2=retcode1; retcode3=retcode1;
end
%how to pick gh and xh
if f3<f & badg3==0
if dispIndx, ih=3, end
fh=f3;xh=x3;gh=g3;badgh=badg3;retcodeh=retcode3;
elseif f2<f & badg2==0
if dispIndx, ih=2, end
fh=f2;xh=x2;gh=g2;badgh=badg2;retcodeh=retcode2;
elseif f1<f & badg1==0
if dispIndx, ih=1, end
fh=f1;xh=x1;gh=g1;badgh=badg1;retcodeh=retcode1;
else
[fh,ih] = min([f1,f2,f3]);
if dispIndx, disp(sprintf('ih = %d',ih)), end
%eval(['xh=x' num2str(ih) ';'])
switch ih
case 1
xh=x1;
case 2
xh=x2;
case 3
xh=x3;
end %case
%eval(['gh=g' num2str(ih) ';'])
%eval(['retcodeh=retcode' num2str(ih) ';'])
retcodei=[retcode1,retcode2,retcode3];
retcodeh=retcodei(ih);
if exist('gh')
nogh=isempty(gh);
else
nogh=1;
end
if nogh
if NumGrad
[gh badgh] = feval('numgrad',fcn,xh,varargin{:});
else
[gh badgh] = feval('grad', xh,varargin{:});
end
end
badgh=1;
end
%end of picking
%ih
%fh
%xh
%gh
%badgh
stuck = (abs(fh-f) < crit);
if (~badg)&(~badgh)&(~stuck)
H = bfgsi(H,gh-g,xh-x);
end
if Verbose
if dispIndx
disp('----')
disp(sprintf('Improvement on iteration %d = %18.9f',itct,f-fh))
end
end
if itct > nit
if dispIndx, disp('iteration count termination'), end
done = 1;
elseif stuck
if dispIndx, disp('improvement < crit termination'), end
done = 1;
end
rc=retcodeh;
if rc == 1
if dispIndx, disp('zero gradient'), end
elseif rc == 6
if dispIndx, disp('smallest step still improving too slow, reversed gradient'), end
elseif rc == 5
if dispIndx, disp('largest step still improving too fast'), end
elseif (rc == 4) | (rc==2)
if dispIndx, disp('back and forth on step length never finished'), end
elseif rc == 3
if dispIndx, disp('smallest step still improving too slow'), end
elseif rc == 7
if dispIndx, disp('warning: possible inaccuracy in H matrix'), end
end
f=fh;
x=xh;
g=gh;
badg=badgh;
end
% what about making an m-file of 10 lines including numgrad.m
% since it appears three times in csminwel.m
matlab/ms-sbvar/cstz/fn_a0freefun.m deleted 100644 → 0
View file @ b00b2fad
function of = fn_a0freefun(b,Ui,nvar,n0,fss,H0inv)
% of = fn_a0freefun(b,Ui,nvar,n0,fss,H0inv)
%
% Negative logPosterior function for squeesed A0 free parameters, which are b's in the WZ notation
% Note: columns correspond to equations
%
% b: sum(n0)-by-1 vector of A0 free parameters
% Ui: nvar-by-1 cell. In each cell, nvar-by-qi orthonormal basis for the null of the ith
% equation contemporaneous restriction matrix where qi is the number of free parameters.
% With this transformation, we have ai = Ui*bi or Ui'*ai = bi where ai is a vector
% of total original parameters and bi is a vector of free parameters. When no
% restrictions are imposed, we have Ui = I. There must be at least one free
% parameter left for the ith equation.
% nvar: number of endogeous variables
% n0: nvar-by-1, ith element represents the number of free A0 parameters in ith equation
% fss: nSample-lags (plus ndobs if dummies are included)
% H0inv: cell(nvar,1). In each cell, posterior inverse of covariance inv(H0) for the ith equation,
% resembling old SpH in the exponent term in posterior of A0, but not divided by T yet.
%----------------
% of: objective function (negative logPosterior)
%
% Tao Zha, February 2000
b=b(:); n0=n0(:);
A0 = zeros(nvar);
n0cum = [0;cumsum(n0)];
tra = 0.0;
for kj = 1:nvar
bj = b(n0cum(kj)+1:n0cum(kj+1));
A0(:,kj) = Ui{kj}*bj;
tra = tra + 0.5*bj'*H0inv{kj}*bj; % negative exponential term
end
[A0l,A0u] = lu(A0);
ada = -fss*sum(log(abs(diag(A0u)))); % negative log determinant of A0 raised to power T
of = ada + tra;
matlab/ms-sbvar/cstz/fn_a0freegrad.m deleted 100644 → 0
View file @ b00b2fad
function [g,badg] = fn_a0freegrad(b,Ui,nvar,n0,fss,H0inv)
% [g,badg] = a0freegrad(b,Ui,nvar,n0,fss,H0inv)
% Analytical gradient for a0freefun.m in use of csminwel.m. See Dhrymes's book.
%
% b: sum(n0)-by-1 vector of A0 free parameters
% Ui: nvar-by-1 cell. In each cell, nvar-by-qi orthonormal basis for the null of the ith
% equation contemporaneous restriction matrix where qi is the number of free parameters.
% With this transformation, we have ai = Ui*bi or Ui'*ai = bi where ai is a vector
% of total original parameters and bi is a vector of free parameters. When no
% restrictions are imposed, we have Ui = I. There must be at least one free
% parameter left for the ith equation.
% nvar: number of endogeous variables
% n0: nvar-by-1, ith element represents the number of free A0 parameters in ith equation
% fss: nSample-lags (plus ndobs if dummies are included)
% H0inv: cell(nvar,1). In each cell, posterior inverse of covariance inv(H0) for the ith equation,
% resembling old SpH in the exponent term in posterior of A0, but not divided by T yet.
%---------------
% g: sum(n0)-by-1 analytical gradient for a0freefun.m
% badg: 0, the value that is used in csminwel.m
%
% Tao Zha, February 2000. Revised, August 2000
b=b(:); n0 = n0(:);
A0 = zeros(nvar);
n0cum = [0;cumsum(n0)];
g = zeros(n0cum(end),1);
badg = 0;
%*** The derivative of the exponential term w.r.t. each free paramater
for kj = 1:nvar
bj = b(n0cum(kj)+1:n0cum(kj+1));
g(n0cum(kj)+1:n0cum(kj+1)) = H0inv{kj}*bj;
A0(:,kj) = Ui{kj}*bj;
end
B=inv(A0');
%*** Add the derivative of -Tlog|A0| w.r.t. each free paramater
for ki = 1:sum(n0)
n = max(find( (ki-n0cum)>0 )); % note, 1<=n<=nvar
g(ki) = g(ki) - fss*B(:,n)'*Ui{n}(:,ki-n0cum(n));
end
matlab/ms-sbvar/cstz/fn_calyrqm.m deleted 100644 → 0
View file @ b00b2fad
function [Myrqm,nMyrqm] = fn_calyrqm(q_m,Byrqm,Eyrqm)
% [Myrqm,nMyrqm] = fn_calyrqm(q_m,Byrqm,Eyrqm)
%
% Given the beginning and end years and quarters (months), export a matrix of all years and
% quarters (months) for these years and in between
%
% q_m: 4 if quarterly and 12 if monthly
% Byrqm: [year quarter(month)] -- all integers, the begining year and quarter (month)
% Eyrqm: [year quarter(month)] -- all integers, the end year and quarter (month)
%-------------------
% Myrqm: matrix of all years and quarters (months) between and incl. Byrqm and Eyrqm
% nMyrqm: number of data points incl. Byrqm and Eyrqm
%
% Tao Zha, April 2000
if ~isempty(find(Byrqm-round(Byrqm))) | (q_m-round(q_m)) | ~isempty(find(Byrqm-round(Byrqm)))
error('argin qm, Byrqm, or Eyrqm must of integer')
elseif Byrqm(1)>Eyrqm(1)
error('Eyrqm(1) must be equal to or greater than Byrqm(1)')
elseif Byrqm(1)==Eyrqm(1)
if Byrqm(2)>Eyrqm(2)
error('Eyrqm(2) must be equal to or greater than Byrqm(2) because of the same year')
end
end
Yr = Byrqm(1)+[0:Eyrqm(1)-Byrqm(1)]';
if length(Yr)>=2 % there are years and quarters (months) between Byrqm and Eyrqm
n=length(Yr)-2;
C=zeros(n*q_m,2);
C(:,1) = kron(Yr(2:end-1),ones(q_m,1));
C(:,2) = kron(ones(n,1),[1:q_m]');
%* initialize a matrix of years and quarters (months) including Byrqm and Eyrqm
Myrqm = zeros((q_m-Byrqm(2)+1)+Eyrqm(2)+n*q_m,2);
%* Years in between
n1=q_m-Byrqm(2)+1;
n2=Eyrqm(2);
Myrqm(n1+1:end-n2,:) = C;
%* Beginning year
for k=1:n1
Myrqm(k,:) = [Byrqm(1) Byrqm(2)+k-1];
end
%* End year
for k=1:n2
Myrqm(end-Eyrqm(2)+k,:) = [Eyrqm(1) k];
end
else %* all the data are in the same calendar year
n1=Eyrqm(2)-Byrqm(2)+1;
Myrqm = zeros(n1,2);
for k=1:n1
Myrqm(k,:) = [Byrqm(1) Byrqm(2)+k-1];
end
end
nMyrqm = size(Myrqm,1);
matlab/ms-sbvar/cstz/fn_dataext.m deleted 100644 → 0
View file @ b00b2fad
function [xdsube,Brow,Erow] = fn_dataext(Byrqm,Eyrqm,xdatae)
% xdsube = dataext(Byrqm,Eyrqm,xdatae)
% Extract subset of xdatae, given Byrqm and Eyrqm
%
% Byrqm: [year quarter(month)]: beginning year and period. If Byqm(2)=0, we get annual data.
% Eyrqm: [year period]: end year and period. If Byrqm(2)=0, it must be Eyrqm(2)=0.
% xdatae: all data (some of which may be NaN) with the first 2 columns indicating years and periods.
%----------
% xdsube: subset of xdatae.
% Brow: the row number in xdatee that marks the first row of xdsube.
% Erow: the row number in xdatee that marks the last row of xdsube.
%
% Tao Zha, April 2000
if (Byrqm(2)==0) & (Eyrqm(2)~=0)
error('If annual data, make sure both Byrqm(2) and Eyrqm(2) are zero')
end
Brow = min(find(xdatae(:,1)==Byrqm(1)));
if isempty(Brow)
error('Byrqm is outside the date range indicated by xdatae(:,1:2)!')
end
if Byrqm(2)>0
nadt=Byrqm(2)-xdatae(Brow,2);
if nadt<0
error('Byrqm is outside the date range indicated by xdatae(:,1:2)!')
end
Brow=Brow+nadt;
end
%
Erow = min(find(xdatae(:,1)==Eyrqm(1)));
if isempty(Erow)
error('Eyrqm is outside the date range indicated by xdatae(:,1:2)!')
end
if Eyrqm(2)>0
nadt2=Eyrqm(2)-xdatae(Erow,2);
if nadt<0
error('Eyrqm is outside the date range indicated by xdatae(:,1:2)!')
end
Erow=Erow+nadt2;
end
%
xdsube = xdatae(Brow:Erow,:); % with dates
matlab/ms-sbvar/cstz/fn_datana.m deleted 100644 → 0
View file @ b00b2fad
function [yactyrge,yactyre,yactqmyge,yactqmge,yactqme] = fn_datana(xdatae,q_m,vlistlog,vlistper,Byrqm,Eyrqm)
% [yactyrge,yactyre,yactqmyge,yactqmge,yactqme] = fn_datana(Byrqm,Eyrqm,xdatae,q_m,vlistlog,vlistper)
%
% Generate prior period, period-to-period-in-last-year, and annual growth rates
% For annual rates, works for both calendar and any annual years, depending on Byrqm and Eyrqm
% If monthly data, we haven't got time to get quarterly growth rates yet.
%
% xdatae: all data (logged levels or interest rates/100, some of which may be NaN) with the first
% 2 columns indicating years and periods.
% q_m: quarter or month period
% vlistlog: sublist for logged variables
% vlistper: sublists for percent variables
% Byrqm: [year quarter(month)]: beginning year and period. If Byqm(2)~=1, we don't get
% calendar annual rates. In other words, the first column of yactyge (which
% indicates years) does not mean calendar years. Byqm(2) must be specified; in other
% words, it must be not set to 0 as in, say, fn_dataext.
% Eyrqm: [year period]: end year and period. Eyqm(2) must be specified; in other words, it
% must be not set to 0 as in, say, fn_dataext.
% NOTE: if no inputs Byrqm and Eyrqm are specified, all growth rates begin at xdatae(1,1:2).
%----------
% yactyrge: annual growth rates with dates in the first 2 columns.
% yactyre: annual average logged level with dates in the 1st 2 columns.
% yactqmyge: period-to-period-in-last-year annual growth rates with dates in the first 2 columns.
% yactqmge: prior-period annualized growth rates with dates in the first 2 columns.
% yactqme: data (logged levels or interest rates/100) with dates in the first 2 columns.
% Same as xdatae but with Brow:Erow.
%
% Tao Zha, April 2000.
if size(xdatae,1)<2*q_m
error('We need at least two years of xdatae to get annual rates. Check xdatae!!')
end
if nargin==4
Brow=1; Erow=length(xdatae(:,1));
nyr = floor((Erow-Brow+1)/q_m);
yrsg = [xdatae(q_m+1,1):xdatae(q_m+1,1)+nyr-2]'; % for annual growth later on
else
if Byrqm(2)<1 | Eyrqm(2)<1
error('This function requires specifying both years and months (quarters) in Byrqm and Eyrqm')
end
Brow = min(find(xdatae(:,1)==Byrqm(1)));
if isempty(Brow)
error('Byrqm is outside the date range of xdatae(:,1:2)!')
end
nadt=Byrqm(2)-xdatae(Brow,2);
if nadt<0
error('Byrqm is outside the date range indicated by xdatae(:,1:2)!')
end
Brow=Brow+nadt;
%
Erow = min(find(xdatae(:,1)==Eyrqm(1)));
if isempty(Brow)
error('Eyrqm is outside the date range of xdatae(:,1:2)!')
end
nadt2=Eyrqm(2)-xdatae(Erow,2);
if nadt<0
error('Eyrqm is outside the date range indicated by xdatae(:,1:2)!')
end
Erow=Erow+nadt2;
nyr = floor((Erow-Brow+1)/q_m);
yrsg = [Byrqm(1)+1:Byrqm(1)+nyr-1]'; % for annual growth later on, which will
% start at Byrqm(1) instead of Byrqm(1)+1
end
%
yactqme = xdatae(Brow:Erow,:); % with dates
yactqm = yactqme(:,3:end); % only data
%======== prior period change (annaluized rate)
yactqmg = yactqm(2:end,:); % start at second period to get growth rate
yactqmg(:,vlistlog) = (yactqm(2:end,vlistlog) - yactqm(1:end-1,vlistlog)) .* q_m;
% monthly, 12*log(1+growth rate), annualized growth rate
yactqmg(:,vlistlog) = 100*(exp(yactqmg(:,vlistlog))-1);
yactqmg(:,vlistper) = 100*yactqmg(:,vlistper);
yactqmge = [yactqme(2:end,1:2) yactqmg];
%======== change from the last year
yactqmyg = yactqm(q_m+1:end,:); % start at the last-year period to get growth rate
yactqmyg(:,vlistlog) = (yactqm(q_m+1:end,vlistlog) - yactqm(1:end-q_m,vlistlog));
yactqmyg(:,vlistlog) = 100*(exp(yactqmyg(:,vlistlog))-1);
yactqmyg(:,vlistper) = 100*yactqmyg(:,vlistper);
yactqmyge = [yactqme(q_m+1:end,1:2) yactqmyg];
%======== annual growth rates
nvar = length(xdatae(1,3:end));
ygmts = yactqm(1:nyr*q_m,:); % converted to the multiplication of q_m
ygmts1 = reshape(ygmts,q_m,nyr,nvar);
ygmts2 = sum(ygmts1,1) ./ q_m;
ygmts3 = reshape(ygmts2,nyr,nvar); % converted to annual average series
%
yactyrg = ygmts3(2:end,:); % start at the last-year period to get growth rate
yactyrg(:,vlistlog) = ygmts3(2:end,vlistlog) - ygmts3(1:end-1,vlistlog);
% annual rate: log(1+growth rate)
yactyrg(:,vlistlog) = 100*(exp(yactyrg(:,vlistlog))-1);
yactyrg(:,vlistper) = 100*yactyrg(:,vlistper);
yactyrge = [yrsg zeros(nyr-1,1) yactyrg];
yrsg1=[yrsg(1)-1:yrsg(end)]';
yactyre = [yrsg1 zeros(nyr,1) ygmts3];
matlab/ms-sbvar/cstz/fn_dataxy.m deleted 100644 → 0
View file @ b00b2fad
function [xtx,xty,yty,fss,phi,y,ncoef,xr,Bh,e] = fn_dataxy(nvar,lags,z,mu,indxDummy,nexo)
% [xtx,xty,yty,fss,phi,y,ncoef,xr,Bh,e] = fn_dataxy(nvar,lags,z,mu,indxDummy,nexo)
%
% Export arranged data matrices for future estimation of the VAR.
% If indxDummy=0, no mu's are used at all and Bh is the OLS estimate.
% If indxDummy=1, only mu(5) and mu(6) as dummy observations. See fn_rnrprior*.m for using mu(1)-mu(4).
% See Wagonner and Zha's Gibbs sampling paper.
%
% nvar: number of endogenous variables.
% lags: the maximum length of lag
% z: T*(nvar+(nexo-1)) matrix of raw or original data (no manipulation involved)
% with sample size including lags and with exogenous variables other than a constant.
% Order of columns: (1) nvar endogenous variables; (2) (nexo-1) exogenous variables;
% (3) constants will be automatically put in the last column.
% mu: 6-by-1 vector of hyperparameters (the following numbers for Atlanta Fed's forecast), where
% only mu(5) and mu(6) are used for dummy observations in this function (i.e.,
% mu(1)-mu(4) are irrelevant here). See fn_rnrprior*.m for using mu(1)-mu(4).
% mu(1): overall tightness and also for A0; (0.57)
% mu(2): relative tightness for A+; (0.13)
% mu(3): relative tightness for the constant term; (0.1)
% mu(4): tightness on lag decay; (1)
% mu(5): weight on nvar sums of coeffs dummy observations (unit roots); (5)
% mu(6): weight on single dummy initial observation including constant
% (cointegration, unit roots, and stationarity); (5)
% indxDummy: 1: add dummy observations to the data; 0: no dummy added.
% nexo: number of exogenous variables. The constant term is the default setting. Besides this term,
% we have nexo-1 exogenous variables. Optional. If left blank, nexo is set to 1.
% -------------------
% xtx: X'X: k-by-k where k=ncoef
% xty: X'Y: k-by-nvar
% yty: Y'Y: nvar-by-nvar
% fss: T: sample size excluding lags. With dummyies, fss=nSample-lags+ndobs.
% phi: X; T-by-k; column: [nvar for 1st lag, ..., nvar for last lag, other exogenous terms, const term]
% y: Y: T-by-nvar where T=fss
% ncoef: number of coefficients in *each* equation. RHS coefficients only, nvar*lags+nexo
% xr: the economy size (ncoef-by-ncoef) in qr(phi) so that xr=chol(X'*X) or xr'*xr=X'*X
% Bh: ncoef-by-nvar estimated reduced-form parameter; column: nvar;
% row: ncoef=[nvar for 1st lag, ..., nvar for last lag, other exogenous terms, const term]
% e: estimated residual e = y -x*Bh, T-by-nvar
%
% Tao Zha, February 2000.
if nargin == 5
nexo=1; % default for constant term
elseif nexo<1
error('We need at least one exogenous term so nexo must >= 1')
end
%*** original sample dimension without dummy prior
nSample = size(z,1); % the sample size (including lags, of course)
sb = lags+1; % original beginning without dummies
ncoef = nvar*lags+nexo; % number of coefficients in *each* equation, RHS coefficients only.
if indxDummy % prior dummy prior
%*** expanded sample dimension by dummy prior
ndobs=nvar+1; % number of dummy observations
fss = nSample+ndobs-lags;
%
% **** nvar prior dummy observations with the sum of coefficients
% ** construct X for Y = X*B + U where phi = X: (T-lags)*k, Y: (T-lags)*nvar
% ** columns: k = # of [nvar for 1st lag, ..., nvar for last lag, exo var, const]
% ** Now, T=T+ndobs -- added with "ndobs" dummy observations
%
phi = zeros(fss,ncoef);
%* constant term
const = ones(fss,1);
const(1:nvar) = zeros(nvar,1);
phi(:,ncoef) = const; % the first nvar periods: no or zero constant!
%* other exogenous (than) constant term
phi(ndobs+1:end,ncoef-nexo+1:ncoef-1) = z(lags+1:end,nvar+1:nvar+nexo-1);
exox = zeros(ndobs,nexo);
phi(1:ndobs,ncoef-nexo+1:ncoef-1) = exox(:,1:nexo-1);
% this = [] when nexo=1 (no other exogenous than constant)
xdgel = z(:,1:nvar); % endogenous variable matrix
xdgelint = mean(xdgel(1:lags,:),1); % mean of the first lags initial conditions
%* Dummies
for k=1:nvar
for m=1:lags
phi(ndobs,nvar*(m-1)+k) = xdgelint(k);
phi(k,nvar*(m-1)+k) = xdgelint(k);
% <<>> multiply hyperparameter later
end
end
%* True data
for k=1:lags
phi(ndobs+1:fss,nvar*(k-1)+1:nvar*k) = xdgel(sb-k:nSample-k,:);
% row: T-lags; column: [nvar for 1st lag, ..., nvar for last lag, exo var, const]
% Thus, # of columns is nvar*lags+nexo = ncoef.
end
%
% ** Y with "ndobs" dummies added
y = zeros(fss,nvar);
%* Dummies
for k=1:nvar
y(ndobs,k) = xdgelint(k);
y(k,k) = xdgelint(k);
% multiply hyperparameter later
end
%* True data
y(ndobs+1:fss,:) = xdgel(sb:nSample,:);
phi(1:nvar,:) = 1*mu(5)*phi(1:nvar,:); % standard Sims and Zha prior
y(1:nvar,:) = mu(5)*y(1:nvar,:); % standard Sims and Zha prior
phi(nvar+1,:) = mu(6)*phi(nvar+1,:);
y(nvar+1,:) = mu(6)*y(nvar+1,:);
[xq,xr]=qr(phi,0);
xtx=xr'*xr;
xty=phi'*y;
[yq,yr]=qr(y,0);
yty=yr'*yr;
Bh = xr\(xr'\xty); % xtx\xty where inv(X'X)*(X'Y)
e=y-phi*Bh;
else
fss = nSample-lags;
%
% ** construct X for Y = X*B + U where phi = X: (T-lags)*k, Y: (T-lags)*nvar
% ** columns: k = # of [nvar for 1st lag, ..., nvar for last lag, exo var, const]
%
phi = zeros(fss,ncoef);
%* constant term
const = ones(fss,1);
phi(:,ncoef) = const; % the first nvar periods: no or zero constant!
%* other exogenous (than) constant term
phi(:,ncoef-nexo+1:ncoef-1) = z(lags+1:end,nvar+1:nvar+nexo-1);
% this = [] when nexo=1 (no other exogenous than constant)
xdgel = z(:,1:nvar); % endogenous variable matrix
%* True data
for k=1:lags
phi(:,nvar*(k-1)+1:nvar*k) = xdgel(sb-k:nSample-k,:);
% row: T-lags; column: [nvar for 1st lag, ..., nvar for last lag, exo var, const]
% Thus, # of columns is nvar*lags+nexo = ncoef.
end
%
y = xdgel(sb:nSample,:);
[xq,xr]=qr(phi,0);
xtx=xr'*xr;
xty=phi'*y;
[yq,yr]=qr(y,0);
yty=yr'*yr;
Bh = xr\(xr'\xty); % xtx\xty where inv(X'X)*(X'Y)
e=y-phi*Bh;
end