-
Notifications
You must be signed in to change notification settings - Fork 6
/
transitionSimBetaBayes.m
72 lines (64 loc) · 2.23 KB
/
transitionSimBetaBayes.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
function [T,C] = transitionSimBetaBayes(nobs, nlag, constant, ...
ks, ys, zs, hs, Vprior, Tprior, Tres)
% transitionSimBayes - Draw parameters from state equation using
% Bayesian methods
%
% Input:
%
% Output:
%
% Usage:
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Copyright (C) 2014 PROFOR Team
%
% 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 3 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/>.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
maxNumberOfRhs = max(ks);
nvary = numel(ks);
stabilityCondition = 1;
while stabilityCondition == 1
% Compute some usefull parameters
xhx = zs'*hs*zs;
xhy = zs'*hs*ys;
capv0inv = inv(Vprior);
%draw from beta conditional on H
capv1inv = capv0inv + xhx;
capv1 = inv(capv1inv);
b1 = capv1*(capv0inv*Tprior + xhy);
beta = b1 + norm_rnd(capv1);
% Transform ALPHA into form that can be used in other utils below and
% stored as output
betan = zeros(nvary, maxNumberOfRhs);
cnt = 0;
for i = 1:nvary
betan(i, Tres(i,:) == 1) = beta(1 + cnt:ks(i) + cnt);
cnt = cnt + ks(i);
end;
% get comp form to check stability
if constant
[T,C] = varGetCompForm(betan(:,2:end,:),betan(:,1,:), nlag, nvary);
else
[T,C] = varGetCompForm(betan, [], nlag, nvary);
end;
if any(abs(eig(T)) >= 1)
stabilityCondition = 1;%1
else
% To brake out of while statement if stability
% ensured
stabilityCondition = 0;
end
end