diff --git a/matlab/ols/sur.m b/matlab/ols/sur.m index c48828944..2d9a466b6 100644 --- a/matlab/ols/sur.m +++ b/matlab/ols/sur.m @@ -11,7 +11,7 @@ function varargout = sur(ds) % SPECIAL REQUIREMENTS % dynare must be run with the option: json=compute -% Copyright (C) 2017 Dynare Team +% Copyright (C) 2017-2018 Dynare Team % % This file is part of Dynare. % @@ -77,10 +77,10 @@ for i = 1:length(lhs) end Y = newY; X = newX; +oo_.sur.dof = length(maxfp:minlp); %% Estimation % Estimated Parameters -oo_.sur.dof = length(maxfp:minlp); [q, r] = qr(X, 0); xpxi = (r'*r)\eye(size(X, 2)); resid = Y - X * (r\(q'*Y)); @@ -90,6 +90,20 @@ M_.Sigma_e = resid'*resid/oo_.sur.dof; kLeye = kron(chol(inv(M_.Sigma_e)), eye(oo_.sur.dof)); [q, r] = qr(kLeye*X, 0); oo_.sur.beta = r\(q'*kLeye*Y); + +%% Return to surgibbs if called from there +st = dbstack(1); +if strcmp(st(1).name, 'surgibbs') + varargout{1} = oo_.sur.dof; + varargout{2} = size(X, 2); + varargout{3} = pidxs; + varargout{4} = oo_.sur.beta; + varargout{5} = X; + varargout{6} = Y; + varargout{7} = length(lhs); + return +end + M_.params(pidxs, 1) = oo_.sur.beta; % Yhat @@ -98,19 +112,6 @@ oo_.sur.Yhat = X * oo_.sur.beta; % Residuals oo_.sur.resid = Y - oo_.sur.Yhat; -%% Return to surgibbs if called from there -st = dbstack(1); -if strcmp(st(1).name, 'surgibbs') - varargout{1} = oo_.sur.dof; - varargout{2} = size(X, 2); - varargout{3} = cellstr(M_.param_names(pidxs, :)); - varargout{4} = oo_.sur.beta; - varargout{5} = X; - varargout{6} = Y; - varargout{7} = length(lhs); - return -end - %% Calculate statistics % Estimate for sigma^2 SS_res = oo_.sur.resid'*oo_.sur.resid; diff --git a/matlab/surgibbs.m b/matlab/surgibbs.m index dbe22a795..e33d56e50 100644 --- a/matlab/surgibbs.m +++ b/matlab/surgibbs.m @@ -1,11 +1,12 @@ -function surgibbs(ds, A, ndraws) -%function surgibbs(ds, A, ndraws) +function surgibbs(ds, A, ndraws, discarddraws) +%function surgibbs(ds, ndraws, discarddraws) % Implements Gibbs Samipling for SUR % % INPUTS -% ds [dseries] data -% A [matrix] prior distribution variance -% ndraws [int] number of draws +% ds [dseries] data +% A [matrix] prior distribution variance +% ndraws [int] number of draws +% discarddraws [int] number of draws to discard % % OUTPUTS % none @@ -13,7 +14,7 @@ function surgibbs(ds, A, ndraws) % SPECIAL REQUIREMENTS % none -% Copyright (C) 2017 Dynare Team +% Copyright (C) 2017-2018 Dynare Team % % This file is part of Dynare. % @@ -30,33 +31,51 @@ function surgibbs(ds, A, ndraws) % You should have received a copy of the GNU General Public License % along with Dynare. If not, see . +%% The notation that follows comes from Section 2.2 of +% Ando, Tomohiro and Zellner, Arnold. 2010. Hierarchical Bayesian Analysis of the +% Seemingly Unrelated Regression and Simultaneous Equations Models Using a +% Combination of Direct Monte Carlo and Importance Sampling Techniques. +% Bayesian Analysis Volume 5, Number 1, pp. 65?96. + global M_ oo_ %% Check input -assert(nargin == 3 || nargin == 5, 'Incorrect number of arguments passed to surgibbs'); +assert(nargin == 3 || nargin == 4, 'Incorrect number of arguments passed to surgibbs'); +if nargin == 3 + discarddraws = 0; +end %% Estimation -[nobs, nvars, pnamesall, beta, X, Y, m] = sur(ds); -beta0 = beta; - -oo_.surgibbs.betadraws = zeros(ndraws, nvars); +beta0 = M_.params; +[nobs, nvars, pidxs, beta, X, Y, m] = sur(ds); +pnamesall = cellstr(M_.param_names(pidxs, :)); +if any(isnan(beta0)) + beta0 = beta; +else + beta = beta0; +end +A = inv(A); +oo_.surgibbs.betadraws = zeros(ndraws-discarddraws, nvars); for i = 1:ndraws - % Draw S, given X, Y, Beta + % Draw Omega, given X, Y, Beta resid = reshape(Y - X*beta, nobs, m); - Omega = rand_inverse_wishart(m, nobs, (resid'*resid)/nobs); + Omega = rand_inverse_wishart(m, nobs, chol(inv(resid'*resid/nobs))); - % Draw beta, given X, Y, S + % Draw beta, given X, Y, Omega tmp = kron(inv(Omega), eye(nobs)); tmp1 = X'*tmp*X; Omegabar = inv(tmp1 + A); - betabar = Omegabar*(tmp1*(tmp1\X'*tmp*Y)+A\beta0); - Sigma_upper_chol = chol(Omegabar, 'upper'); - beta = rand_multivariate_normal(betabar', Sigma_upper_chol, nvars)'; - oo_.surgibbs.betadraws(i, :) = beta'; + betahat = tmp1\X'*tmp*Y; + betabar = Omegabar*(tmp1*betahat+A*beta0); + beta = rand_multivariate_normal(betabar', chol(Omegabar), nvars)'; + if i > discarddraws + oo_.surgibbs.betadraws(i-discarddraws, :) = beta'; + end end % save parameter values -oo_.surgibbs.beta = (sum(oo_.surgibbs.betadraws)/ndraws)'; +oo_.surgibbs.beta = (sum(oo_.surgibbs.betadraws)/(ndraws-discarddraws))'; +M_.params(pidxs, 1) = oo_.surgibbs.beta; %% Print Output dyn_table('Gibbs Sampling on SUR', {}, {}, pnamesall, ... @@ -77,4 +96,3 @@ for j = 1:length(pnamesall) line([oo_.surgibbs.beta(j) oo_.surgibbs.beta(j)], [min(hc) max(hc)], 'Color', 'red'); title(pnamesall{j}, 'Interpreter', 'none') end - diff --git a/tests/ECB/SURGibbs/fishdata.csv b/tests/ECB/SURGibbs/fishdata.csv new file mode 100644 index 000000000..5d525ab7a --- /dev/null +++ b/tests/ECB/SURGibbs/fishdata.csv @@ -0,0 +1,112 @@ +day1,day2,day3,day4,date,stormy,mixed,price,qty,rainy,cold,windspd,windspd2,pricelevel,totr,tots +1,0,0,0,911202,1,0,-0.4307829,8.994421,1,0,2.995732,8.974412,0.65,7232,8058 +0,1,0,0,911203,1,0,0,7.707063,0,0,2.995732,8.974412,1,2110,2224 +0,0,1,0,911204,0,1,0.0723207,8.350194,1,1,2.813411,7.91528,1.075,5247,4231 +0,0,0,1,911205,1,0,0.247139,8.656955,0,1,3.036554,9.220662,1.280357,1290,5750 +0,0,0,0,911206,1,0,0.6643268,7.844241,0,1,3.036554,9.220662,1.943182,1717,2551 +1,0,0,0,911209,0,0,-0.2065143,9.301277,0,0,2.762117,7.629292,0.8134146,11643,10977 +0,1,0,0,911210,0,1,-0.1158318,8.920656,0,0,2.762117,7.629292,0.890625,9640,7485 +0,0,1,0,911211,0,0,-0.2598674,9.105979,1,0,2.762117,7.629292,0.7711539,9347,9009 +0,0,0,1,911212,0,1,-0.1171254,8.307706,0,0,2.762117,7.629292,0.8894737,3890,4055 +0,0,0,0,911213,0,0,-0.3420761,9.20954,0,0,2.590267,6.709484,0.7102941,16318,9992 +1,0,0,0,911216,1,0,-0.1255632,8.552561,0,1,3.149883,9.921763,0.882,8725,5180 +0,1,0,0,911217,1,0,0.027399,8.523175,0,1,3.218876,10.36116,1.027778,2780,5030 +0,0,1,0,911218,1,0,-0.0712275,8.865453,0,1,3.184974,10.14406,0.93125,9078,6833 +0,0,0,1,911219,1,0,0.1230601,9.186355,0,1,3.113515,9.693978,1.130952,5066,9763 +0,0,0,0,911220,1,0,0.2130932,8.699348,0,1,3.149883,9.921763,1.2375,4796,5999 +1,0,0,0,911223,0,1,-0.3172045,9.408863,0,1,2.953173,8.721229,0.7281818,13647,12196 +0,1,0,0,911224,0,1,-0.1088388,8.150179,0,1,3.036554,9.220662,0.896875,1255,3464 +0,0,0,1,911226,0,1,0.2231435,6.703188,0,1,2.953173,8.721229,1.25,1115,815 +0,0,0,0,911227,0,0,0.2464593,8.798907,0,1,2.862201,8.192194,1.279487,6887,6627 +1,0,0,0,911230,0,0,-0.075431,9.565214,1,1,2.70805,7.333536,0.9273437,15894,14260 +0,1,0,0,911231,0,0,0.2055992,8.297792,1,1,2.862201,8.192194,1.228261,5850,4015 +0,0,0,1,920102,0,0,0.2188098,8.320935,0,0,2.813411,7.91528,1.244595,409,4109 +0,0,0,0,920103,0,0,0.307025,8.884887,0,0,2.456736,6.035551,1.359375,7222,7222 +1,0,0,0,920106,1,0,0.399592,9.336444,0,0,3.314186,10.98383,1.491216,13036,11224 +0,1,0,0,920107,1,0,0.4969802,8.122668,0,1,3.344039,11.1826,1.64375,1760,3370 +0,0,1,0,920108,1,0,0.3968258,8.15191,0,1,3.036554,9.220662,1.487097,4824,3470 +0,0,0,1,920109,0,1,0.2830518,9.51834,0,1,2.813411,7.91528,1.327174,16489,13607 +0,0,0,0,920110,0,1,0.2263384,8.567886,1,0,2.650892,7.027227,1.254,4842,5260 +1,0,0,0,920113,0,1,-0.0723207,9.386811,0,0,2.908721,8.460658,0.9302325,12732,11930 +0,1,0,0,920114,0,1,-0.1545295,8.628735,0,0,2.813411,7.91528,0.8568182,7070,5590 +0,0,1,0,920115,1,0,0.2363888,7.727535,1,1,2.995732,8.974412,1.266667,2873,2270 +0,0,0,1,920116,1,0,0.1718503,8.453188,0,1,3.113515,9.693978,1.1875,1915,4690 +0,0,0,0,920117,1,0,0.1541506,7.156956,1,1,3.218876,10.36116,1.166667,240,1283 +1,0,0,0,920120,1,0,0.0350913,6.363028,0,1,3.075775,9.460391,1.035714,300,580 +0,1,0,0,920121,0,1,0.218131,7.313221,1,1,3.075775,9.460391,1.24375,1960,1550 +0,0,1,0,920122,0,0,-0.0624355,8.39163,0,1,2.908721,8.460658,0.9394737,5408,4410 +0,0,0,1,920123,0,0,-0.4418328,9.188503,0,0,2.70805,7.333536,0.6428571,10130,9909 +0,0,0,0,920124,0,0,-0.7699343,9.318297,1,0,2.862201,8.192194,0.4630435,12943,11140 +1,0,0,0,920127,1,0,-0.7391911,8.363809,1,1,3.113515,9.693978,0.4775,2766,4289 +0,1,0,0,920128,0,1,-0.7230001,8.621193,0,1,2.862201,8.192194,0.4852941,4050,5548 +0,0,1,0,920129,0,0,-0.5184302,8.144679,0,1,2.813411,7.91528,0.5954546,6208,3445 +0,0,0,1,920130,0,0,-0.5533853,8.710784,0,1,2.650892,7.027227,0.575,7539,6151 +0,0,0,0,920131,0,0,-0.5590724,8.417373,0,0,2.650892,7.027227,0.5717391,6250,4525 +1,0,0,0,920203,1,0,-0.4238143,9.017968,0,1,3.251666,10.57333,0.6545454,8620,8250 +0,1,0,0,920204,1,0,-0.7731899,7.528332,0,1,3.283414,10.78081,0.4615385,300,1860 +0,0,1,0,920205,1,0,-0.2040953,7.154615,0,1,3.075775,9.460391,0.8153846,360,1280 +0,0,0,1,920206,1,0,0.3841877,8.487764,0,1,3.075775,9.460391,1.468421,5680,4855 +0,0,0,0,920207,0,1,0.3479887,8.979543,0,1,2.953173,8.721229,1.416216,7484,8089 +1,0,0,0,920210,1,0,0.3229889,8.744328,0,1,3.036554,9.220662,1.38125,6960,6275 +0,1,0,0,920211,1,0,0.4345259,7.893572,0,1,2.953173,8.721229,1.544231,2460,2680 +0,0,1,0,920212,1,0,0.2000211,8.5923,0,1,2.908721,8.460658,1.221429,8396,5390 +0,0,0,1,920213,0,1,-0.0071685,9.035987,0,1,2.862201,8.192194,0.9928572,14440,8400 +0,0,0,0,920214,0,1,-0.4744579,9.525735,0,1,2.862201,8.192194,0.6222222,7095,13708 +0,1,0,0,920218,0,0,-0.7711087,9.773094,0,1,2.862201,8.192194,0.4625,19698,17555 +0,0,1,0,920219,0,1,-1.094774,8.158802,0,0,2.813411,7.91528,0.3346154,10,3474 +0,0,0,1,920220,0,0,-0.8602012,8.677609,0,0,2.70805,7.333536,0.4230769,6470,5870 +0,0,0,0,920221,0,1,-0.597837,9.034796,0,0,2.762117,7.629292,0.55,6695,8390 +1,0,0,0,920224,0,0,-0.567984,8.375629,0,1,2.525729,6.379305,0.5666667,5222,4340 +0,1,0,0,920225,0,0,-0.7711087,7.649693,0,1,2.590267,6.709484,0.4625,2885,2100 +0,0,1,0,920226,0,1,-0.2724146,7.853605,0,1,2.862201,8.192194,0.7615384,2025,2575 +0,0,0,1,920227,1,0,0.1849223,7.955074,1,1,2.953173,8.721229,1.203125,1179,2850 +0,0,0,0,920228,1,0,0.20676,8.73069,0,0,2.908721,8.460658,1.229687,7420,6190 +1,0,0,0,920302,0,1,0.074108,8.633197,0,0,2.908721,8.460658,1.076923,5855,5635 +0,1,0,0,920303,1,0,-0.0071685,8.796339,0,1,2.908721,8.460658,0.9928572,7327,6730 +0,0,1,0,920304,0,1,-0.2548923,8.994048,0,0,2.762117,7.629292,0.775,7850,8055 +0,0,0,1,920305,0,1,-0.2364898,9.22375,0,0,2.650892,7.027227,0.789394,8610,10135 +0,0,0,0,920306,0,0,-0.4848501,8.568457,0,0,2.590267,6.709484,0.6157895,5213,5263 +1,0,0,0,920309,0,1,-0.8119307,9.981374,0,0,2.862201,8.192194,0.444,26530,21620 +0,1,0,0,920310,0,1,-0.4780358,9.172119,0,0,2.813411,7.91528,0.62,4520,9685 +0,0,1,0,920311,0,1,-0.2301122,8.116715,0,1,2.762117,7.629292,0.7944444,3375,3350 +0,0,0,1,920312,1,0,-0.2887938,9.241354,1,1,3.075775,9.460391,0.7491667,10600,11015 +0,0,0,0,920313,1,0,-0.3409266,8.016318,0,1,3.283414,10.78081,0.7111111,21070,3030 +1,0,0,0,920316,0,0,-0.4812668,9.646011,0,1,2.862201,8.192194,0.618,250,15520 +0,1,0,0,920317,0,0,-0.5817336,9.168476,0,1,2.813411,7.91528,0.5589285,7610,9590 +0,0,1,0,920318,0,1,-0.0370413,7.44132,0,1,2.908721,8.460658,0.9636363,3065,1705 +0,0,0,1,920319,0,0,0.0746435,8.61432,0,1,2.813411,7.91528,1.0775,3830,5510 +0,0,0,0,920320,1,0,-0.5268259,8.902047,1,1,3.075775,9.460391,0.5904762,8342,7347 +1,0,0,0,920323,0,0,-0.1123292,7.857481,0,1,2.590267,6.709484,0.89375,2305,2585 +0,1,0,0,920324,0,1,-0.2376716,8.228177,0,1,2.762117,7.629292,0.7884616,4410,3745 +0,0,1,0,920325,0,1,-0.4243931,7.986165,0,1,2.762117,7.629292,0.6541666,3000,2940 +0,0,0,1,920326,0,1,-0.4828126,9.395989,0,0,2.862201,8.192194,0.6170455,11783,12140 +0,0,0,0,920327,0,1,-0.1923718,7.711997,0,0,2.862201,8.192194,0.825,2535,2295 +1,0,0,0,920330,1,0,0.1823216,7.705263,0,0,3.401197,11.56814,1.2,5130,2220 +0,1,0,0,920331,1,0,0.2876821,7.013916,0,0,3.184974,10.14406,1.333333,514,1112 +0,0,1,0,920401,0,1,-0.1873341,8.492901,1,0,2.953173,8.721229,0.8291667,3920,4940 +0,0,0,1,920402,0,0,-0.4404448,9.333443,0,0,2.590267,6.709484,0.64375,10240,11310 +0,0,0,0,920403,0,0,-0.567984,9.578795,1,1,2.862201,8.192194,0.5666667,15042,14480 +1,0,0,0,920406,0,0,-0.3873608,8.975504,0,0,2.813411,7.91528,0.6788461,6852,7727 +0,1,0,0,920407,0,0,-0.5108256,7.28961,0,0,2.70805,7.333536,0.6,305,1465 +0,0,1,0,920408,0,0,-0.6149098,9.415728,0,0,2.70805,7.333536,0.5406896,14148,12280 +0,0,0,1,920409,0,0,-0.3876929,8.740336,0,0,2.650892,7.027227,0.6786207,5160,6250 +0,0,0,0,920410,0,0,-0.3716677,8.584852,0,0,2.650892,7.027227,0.6895834,4340,5350 +1,0,0,0,920413,0,0,-0.3101549,8.657825,0,1,2.862201,8.192194,0.7333333,7049,5535 +0,1,0,0,920414,0,0,-0.3502024,7.383989,0,0,2.995732,8.974412,0.7045454,200,2260 +0,0,1,0,920415,0,0,-0.2231435,6.194406,0,0,2.762117,7.629292,0.8,360,490 +0,0,0,1,920416,0,0,-0.3119048,9.187583,0,0,2.525729,6.379305,0.7320513,11875,9655 +0,0,0,0,920417,0,0,-0.2258646,8.779557,1,0,2.456736,6.035551,0.7978261,3900,6050 +0,1,0,0,920421,0,0,-0.2842691,9.261033,0,0,2.590267,6.709484,0.7525641,10700,10494 +0,0,1,0,920422,0,0,-0.0971638,8.294049,0,0,2.525729,6.379305,0.9074074,3940,4118 +0,0,0,1,920423,0,0,-0.1419181,8.737132,1,0,2.650892,7.027227,0.8676924,6570,5910 +0,0,0,0,920424,0,0,-0.5364382,9.651816,0,0,2.650892,7.027227,0.5848276,15940,15455 +1,0,0,0,920427,0,0,-0.4596867,8.552561,0,0,2.650892,7.027227,0.6314815,8170,5060 +0,1,0,0,920428,0,0,-0.4866114,8.163371,0,0,2.590267,6.709484,0.6147059,1676,3300 +0,0,1,0,920429,0,0,-0.7603559,8.663024,0,0,2.525729,6.379305,0.4675,6963,5785 +0,0,0,1,920430,0,0,-0.8931236,9.099409,0,0,2.525729,6.379305,0.409375,7694,9050 +0,0,0,0,920501,0,0,-1.107745,9.014325,0,0,2.525729,6.379305,0.330303,11370,8520 +1,0,0,0,920504,0,0,-0.7985078,8.610683,0,0,2.862201,8.192194,0.45,1610,6260 +0,1,0,0,920505,0,1,-0.0870114,7.162397,0,0,2.908721,8.460658,0.9166666,1320,1290 +0,0,1,0,920506,0,1,0.1849223,7.36201,0,0,2.862201,8.192194,1.203125,2305,1450 +0,0,0,1,920507,0,1,0.2231435,8.764053,0,0,2.813411,7.91528,1.25,6080,4780 +0,0,0,0,920508,0,1,0.5611184,8.328451,0,0,2.862201,8.192194,1.752632,4020,3420 diff --git a/tests/ECB/SURGibbs/fulton_fish.mod b/tests/ECB/SURGibbs/fulton_fish.mod new file mode 100644 index 000000000..f0c4847b2 --- /dev/null +++ b/tests/ECB/SURGibbs/fulton_fish.mod @@ -0,0 +1,24 @@ +// --+ options: json=compute +-- +% Example from Section 6.1 of +% Ando, Tomohiro and Zellner, Arnold. 2010. Hierarchical Bayesian Analysis of the +% Seemingly Unrelated Regression and Simultaneous Equations Models Using a +% Combination of Direct Monte Carlo and Importance Sampling Techniques. +% Bayesian Analysis Volume 5, Number 1, pp. 65?96. + +var qty, price; +varexo res_u, res_v, stormy, mixed; +parameters bq0, bp0, bq1, bp1, bp2; + +bp0 = 8.5527; +bp1 = -0.5302; +bp2 = -0.3974; + +bq0 = 6.7523; +bq1 = -0.7969; + +model(linear); + price = bp0 + bp1*stormy + bp2*mixed + res_v; + qty = bq0 + bq1*price + res_u; +end; + +surgibbs(dseries('fishdata.csv'), 0.0005.*eye(M_.param_nbr), 10000, 5000);