surgibbs: draft change and test

time-shift
Houtan 2018-01-05 17:13:46 +01:00
parent 27819954d0
commit c1de50bbc2
4 changed files with 190 additions and 35 deletions

View File

@ -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;

View File

@ -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 <http://www.gnu.org/licenses/>.
%% 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

View File

@ -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
1 day1 day2 day3 day4 date stormy mixed price qty rainy cold windspd windspd2 pricelevel totr tots
2 1 0 0 0 911202 1 0 -0.4307829 8.994421 1 0 2.995732 8.974412 0.65 7232 8058
3 0 1 0 0 911203 1 0 0 7.707063 0 0 2.995732 8.974412 1 2110 2224
4 0 0 1 0 911204 0 1 0.0723207 8.350194 1 1 2.813411 7.91528 1.075 5247 4231
5 0 0 0 1 911205 1 0 0.247139 8.656955 0 1 3.036554 9.220662 1.280357 1290 5750
6 0 0 0 0 911206 1 0 0.6643268 7.844241 0 1 3.036554 9.220662 1.943182 1717 2551
7 1 0 0 0 911209 0 0 -0.2065143 9.301277 0 0 2.762117 7.629292 0.8134146 11643 10977
8 0 1 0 0 911210 0 1 -0.1158318 8.920656 0 0 2.762117 7.629292 0.890625 9640 7485
9 0 0 1 0 911211 0 0 -0.2598674 9.105979 1 0 2.762117 7.629292 0.7711539 9347 9009
10 0 0 0 1 911212 0 1 -0.1171254 8.307706 0 0 2.762117 7.629292 0.8894737 3890 4055
11 0 0 0 0 911213 0 0 -0.3420761 9.20954 0 0 2.590267 6.709484 0.7102941 16318 9992
12 1 0 0 0 911216 1 0 -0.1255632 8.552561 0 1 3.149883 9.921763 0.882 8725 5180
13 0 1 0 0 911217 1 0 0.027399 8.523175 0 1 3.218876 10.36116 1.027778 2780 5030
14 0 0 1 0 911218 1 0 -0.0712275 8.865453 0 1 3.184974 10.14406 0.93125 9078 6833
15 0 0 0 1 911219 1 0 0.1230601 9.186355 0 1 3.113515 9.693978 1.130952 5066 9763
16 0 0 0 0 911220 1 0 0.2130932 8.699348 0 1 3.149883 9.921763 1.2375 4796 5999
17 1 0 0 0 911223 0 1 -0.3172045 9.408863 0 1 2.953173 8.721229 0.7281818 13647 12196
18 0 1 0 0 911224 0 1 -0.1088388 8.150179 0 1 3.036554 9.220662 0.896875 1255 3464
19 0 0 0 1 911226 0 1 0.2231435 6.703188 0 1 2.953173 8.721229 1.25 1115 815
20 0 0 0 0 911227 0 0 0.2464593 8.798907 0 1 2.862201 8.192194 1.279487 6887 6627
21 1 0 0 0 911230 0 0 -0.075431 9.565214 1 1 2.70805 7.333536 0.9273437 15894 14260
22 0 1 0 0 911231 0 0 0.2055992 8.297792 1 1 2.862201 8.192194 1.228261 5850 4015
23 0 0 0 1 920102 0 0 0.2188098 8.320935 0 0 2.813411 7.91528 1.244595 409 4109
24 0 0 0 0 920103 0 0 0.307025 8.884887 0 0 2.456736 6.035551 1.359375 7222 7222
25 1 0 0 0 920106 1 0 0.399592 9.336444 0 0 3.314186 10.98383 1.491216 13036 11224
26 0 1 0 0 920107 1 0 0.4969802 8.122668 0 1 3.344039 11.1826 1.64375 1760 3370
27 0 0 1 0 920108 1 0 0.3968258 8.15191 0 1 3.036554 9.220662 1.487097 4824 3470
28 0 0 0 1 920109 0 1 0.2830518 9.51834 0 1 2.813411 7.91528 1.327174 16489 13607
29 0 0 0 0 920110 0 1 0.2263384 8.567886 1 0 2.650892 7.027227 1.254 4842 5260
30 1 0 0 0 920113 0 1 -0.0723207 9.386811 0 0 2.908721 8.460658 0.9302325 12732 11930
31 0 1 0 0 920114 0 1 -0.1545295 8.628735 0 0 2.813411 7.91528 0.8568182 7070 5590
32 0 0 1 0 920115 1 0 0.2363888 7.727535 1 1 2.995732 8.974412 1.266667 2873 2270
33 0 0 0 1 920116 1 0 0.1718503 8.453188 0 1 3.113515 9.693978 1.1875 1915 4690
34 0 0 0 0 920117 1 0 0.1541506 7.156956 1 1 3.218876 10.36116 1.166667 240 1283
35 1 0 0 0 920120 1 0 0.0350913 6.363028 0 1 3.075775 9.460391 1.035714 300 580
36 0 1 0 0 920121 0 1 0.218131 7.313221 1 1 3.075775 9.460391 1.24375 1960 1550
37 0 0 1 0 920122 0 0 -0.0624355 8.39163 0 1 2.908721 8.460658 0.9394737 5408 4410
38 0 0 0 1 920123 0 0 -0.4418328 9.188503 0 0 2.70805 7.333536 0.6428571 10130 9909
39 0 0 0 0 920124 0 0 -0.7699343 9.318297 1 0 2.862201 8.192194 0.4630435 12943 11140
40 1 0 0 0 920127 1 0 -0.7391911 8.363809 1 1 3.113515 9.693978 0.4775 2766 4289
41 0 1 0 0 920128 0 1 -0.7230001 8.621193 0 1 2.862201 8.192194 0.4852941 4050 5548
42 0 0 1 0 920129 0 0 -0.5184302 8.144679 0 1 2.813411 7.91528 0.5954546 6208 3445
43 0 0 0 1 920130 0 0 -0.5533853 8.710784 0 1 2.650892 7.027227 0.575 7539 6151
44 0 0 0 0 920131 0 0 -0.5590724 8.417373 0 0 2.650892 7.027227 0.5717391 6250 4525
45 1 0 0 0 920203 1 0 -0.4238143 9.017968 0 1 3.251666 10.57333 0.6545454 8620 8250
46 0 1 0 0 920204 1 0 -0.7731899 7.528332 0 1 3.283414 10.78081 0.4615385 300 1860
47 0 0 1 0 920205 1 0 -0.2040953 7.154615 0 1 3.075775 9.460391 0.8153846 360 1280
48 0 0 0 1 920206 1 0 0.3841877 8.487764 0 1 3.075775 9.460391 1.468421 5680 4855
49 0 0 0 0 920207 0 1 0.3479887 8.979543 0 1 2.953173 8.721229 1.416216 7484 8089
50 1 0 0 0 920210 1 0 0.3229889 8.744328 0 1 3.036554 9.220662 1.38125 6960 6275
51 0 1 0 0 920211 1 0 0.4345259 7.893572 0 1 2.953173 8.721229 1.544231 2460 2680
52 0 0 1 0 920212 1 0 0.2000211 8.5923 0 1 2.908721 8.460658 1.221429 8396 5390
53 0 0 0 1 920213 0 1 -0.0071685 9.035987 0 1 2.862201 8.192194 0.9928572 14440 8400
54 0 0 0 0 920214 0 1 -0.4744579 9.525735 0 1 2.862201 8.192194 0.6222222 7095 13708
55 0 1 0 0 920218 0 0 -0.7711087 9.773094 0 1 2.862201 8.192194 0.4625 19698 17555
56 0 0 1 0 920219 0 1 -1.094774 8.158802 0 0 2.813411 7.91528 0.3346154 10 3474
57 0 0 0 1 920220 0 0 -0.8602012 8.677609 0 0 2.70805 7.333536 0.4230769 6470 5870
58 0 0 0 0 920221 0 1 -0.597837 9.034796 0 0 2.762117 7.629292 0.55 6695 8390
59 1 0 0 0 920224 0 0 -0.567984 8.375629 0 1 2.525729 6.379305 0.5666667 5222 4340
60 0 1 0 0 920225 0 0 -0.7711087 7.649693 0 1 2.590267 6.709484 0.4625 2885 2100
61 0 0 1 0 920226 0 1 -0.2724146 7.853605 0 1 2.862201 8.192194 0.7615384 2025 2575
62 0 0 0 1 920227 1 0 0.1849223 7.955074 1 1 2.953173 8.721229 1.203125 1179 2850
63 0 0 0 0 920228 1 0 0.20676 8.73069 0 0 2.908721 8.460658 1.229687 7420 6190
64 1 0 0 0 920302 0 1 0.074108 8.633197 0 0 2.908721 8.460658 1.076923 5855 5635
65 0 1 0 0 920303 1 0 -0.0071685 8.796339 0 1 2.908721 8.460658 0.9928572 7327 6730
66 0 0 1 0 920304 0 1 -0.2548923 8.994048 0 0 2.762117 7.629292 0.775 7850 8055
67 0 0 0 1 920305 0 1 -0.2364898 9.22375 0 0 2.650892 7.027227 0.789394 8610 10135
68 0 0 0 0 920306 0 0 -0.4848501 8.568457 0 0 2.590267 6.709484 0.6157895 5213 5263
69 1 0 0 0 920309 0 1 -0.8119307 9.981374 0 0 2.862201 8.192194 0.444 26530 21620
70 0 1 0 0 920310 0 1 -0.4780358 9.172119 0 0 2.813411 7.91528 0.62 4520 9685
71 0 0 1 0 920311 0 1 -0.2301122 8.116715 0 1 2.762117 7.629292 0.7944444 3375 3350
72 0 0 0 1 920312 1 0 -0.2887938 9.241354 1 1 3.075775 9.460391 0.7491667 10600 11015
73 0 0 0 0 920313 1 0 -0.3409266 8.016318 0 1 3.283414 10.78081 0.7111111 21070 3030
74 1 0 0 0 920316 0 0 -0.4812668 9.646011 0 1 2.862201 8.192194 0.618 250 15520
75 0 1 0 0 920317 0 0 -0.5817336 9.168476 0 1 2.813411 7.91528 0.5589285 7610 9590
76 0 0 1 0 920318 0 1 -0.0370413 7.44132 0 1 2.908721 8.460658 0.9636363 3065 1705
77 0 0 0 1 920319 0 0 0.0746435 8.61432 0 1 2.813411 7.91528 1.0775 3830 5510
78 0 0 0 0 920320 1 0 -0.5268259 8.902047 1 1 3.075775 9.460391 0.5904762 8342 7347
79 1 0 0 0 920323 0 0 -0.1123292 7.857481 0 1 2.590267 6.709484 0.89375 2305 2585
80 0 1 0 0 920324 0 1 -0.2376716 8.228177 0 1 2.762117 7.629292 0.7884616 4410 3745
81 0 0 1 0 920325 0 1 -0.4243931 7.986165 0 1 2.762117 7.629292 0.6541666 3000 2940
82 0 0 0 1 920326 0 1 -0.4828126 9.395989 0 0 2.862201 8.192194 0.6170455 11783 12140
83 0 0 0 0 920327 0 1 -0.1923718 7.711997 0 0 2.862201 8.192194 0.825 2535 2295
84 1 0 0 0 920330 1 0 0.1823216 7.705263 0 0 3.401197 11.56814 1.2 5130 2220
85 0 1 0 0 920331 1 0 0.2876821 7.013916 0 0 3.184974 10.14406 1.333333 514 1112
86 0 0 1 0 920401 0 1 -0.1873341 8.492901 1 0 2.953173 8.721229 0.8291667 3920 4940
87 0 0 0 1 920402 0 0 -0.4404448 9.333443 0 0 2.590267 6.709484 0.64375 10240 11310
88 0 0 0 0 920403 0 0 -0.567984 9.578795 1 1 2.862201 8.192194 0.5666667 15042 14480
89 1 0 0 0 920406 0 0 -0.3873608 8.975504 0 0 2.813411 7.91528 0.6788461 6852 7727
90 0 1 0 0 920407 0 0 -0.5108256 7.28961 0 0 2.70805 7.333536 0.6 305 1465
91 0 0 1 0 920408 0 0 -0.6149098 9.415728 0 0 2.70805 7.333536 0.5406896 14148 12280
92 0 0 0 1 920409 0 0 -0.3876929 8.740336 0 0 2.650892 7.027227 0.6786207 5160 6250
93 0 0 0 0 920410 0 0 -0.3716677 8.584852 0 0 2.650892 7.027227 0.6895834 4340 5350
94 1 0 0 0 920413 0 0 -0.3101549 8.657825 0 1 2.862201 8.192194 0.7333333 7049 5535
95 0 1 0 0 920414 0 0 -0.3502024 7.383989 0 0 2.995732 8.974412 0.7045454 200 2260
96 0 0 1 0 920415 0 0 -0.2231435 6.194406 0 0 2.762117 7.629292 0.8 360 490
97 0 0 0 1 920416 0 0 -0.3119048 9.187583 0 0 2.525729 6.379305 0.7320513 11875 9655
98 0 0 0 0 920417 0 0 -0.2258646 8.779557 1 0 2.456736 6.035551 0.7978261 3900 6050
99 0 1 0 0 920421 0 0 -0.2842691 9.261033 0 0 2.590267 6.709484 0.7525641 10700 10494
100 0 0 1 0 920422 0 0 -0.0971638 8.294049 0 0 2.525729 6.379305 0.9074074 3940 4118
101 0 0 0 1 920423 0 0 -0.1419181 8.737132 1 0 2.650892 7.027227 0.8676924 6570 5910
102 0 0 0 0 920424 0 0 -0.5364382 9.651816 0 0 2.650892 7.027227 0.5848276 15940 15455
103 1 0 0 0 920427 0 0 -0.4596867 8.552561 0 0 2.650892 7.027227 0.6314815 8170 5060
104 0 1 0 0 920428 0 0 -0.4866114 8.163371 0 0 2.590267 6.709484 0.6147059 1676 3300
105 0 0 1 0 920429 0 0 -0.7603559 8.663024 0 0 2.525729 6.379305 0.4675 6963 5785
106 0 0 0 1 920430 0 0 -0.8931236 9.099409 0 0 2.525729 6.379305 0.409375 7694 9050
107 0 0 0 0 920501 0 0 -1.107745 9.014325 0 0 2.525729 6.379305 0.330303 11370 8520
108 1 0 0 0 920504 0 0 -0.7985078 8.610683 0 0 2.862201 8.192194 0.45 1610 6260
109 0 1 0 0 920505 0 1 -0.0870114 7.162397 0 0 2.908721 8.460658 0.9166666 1320 1290
110 0 0 1 0 920506 0 1 0.1849223 7.36201 0 0 2.862201 8.192194 1.203125 2305 1450
111 0 0 0 1 920507 0 1 0.2231435 8.764053 0 0 2.813411 7.91528 1.25 6080 4780
112 0 0 0 0 920508 0 1 0.5611184 8.328451 0 0 2.862201 8.192194 1.752632 4020 3420

View File

@ -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);