forked from burakbayramli/books
-
Notifications
You must be signed in to change notification settings - Fork 0
/
ch16q3.m
212 lines (185 loc) · 5.24 KB
/
ch16q3.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
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
%Program for Exercise 3 in Chapter 16
%Bayesian model averaging program
%Use cross country growth example from Fernandez, Ley and Steel (JAE, 2001)
load growth.dat;
%The data set is arranged with data for each country taking up 6 lines
%The following makes it into an N by K matrix
n=72;
%There probably is a better way of doing this...
rawdat=zeros(n,42);
j=1;
for i=1:n
rawdat(i,:)= [growth(j,:) growth(j+1,:) growth(j+2,:) ...
growth(j+3,:) growth(j+4,:) growth(j+5,:)];
j=j+6;
end
y=rawdat(:,1);
xraw=rawdat(:,2:42);
bigk=size(xraw,2);
%subtract mean from all regressors as in FLS
mxraw=mean(xraw);
for i=1:bigk
xraw(:,i)=xraw(:,i) - mxraw(1,i);
end
%I am creating a vector of length bigk which is 1 if explanatory
%variable is included, else equals zero
%molddraw is initialized to contain all explanatory variables
%with t-stats greater than .5
molddraw=zeros(bigk,1);
xold=[ones(n,1) xraw];
xtxinv=inv(xold'*xold);
bhat=xtxinv*xold'*y;
e=y-xold*bhat;
sse=e'*e;
s2=sse/(n-bigk-1);
bcov=s2*xtxinv;
bt=zeros(bigk+1,1);
for i=2:bigk+1
bt(i,1)=bhat(i,1)/sqrt(bcov(i,i));
if abs(bt(i,1))>.5
molddraw(i-1,1)=1;
end
end
%Make up the matrix of explanatory variables for model
xold=ones(n,1);
kold=sum(molddraw)+1;
for i=1:bigk
if molddraw(i,1)>0
xold=[xold xraw(:,i)];
end
end
%specify g0 for the g-prior
if n<=(bigk^2)
g0=1/(bigk^2);
else
g0=1/n;
end
yty = (y-mean(y))'*(y-mean(y));
xtxinv=inv(xold'*xold);
ymy = y'*y - y'*xold*xtxinv*xold'*y;
g1=g0/(g0+1);
g2=1/(g0+1);
lprobold = .5*kold*log(g1) -.5*(n-1)*log(g2*ymy + g1*yty);
mstart=molddraw;
lprobstart=lprobold;
inccount=zeros(bigk,1);
msize=0;
%I am keeping records for top 10 drawn models
%Here initialize this to initial model
top10mod=[molddraw molddraw molddraw molddraw molddraw ...
molddraw molddraw molddraw molddraw molddraw];
lprobtop10=lprobold*ones(10,1);
top10count=zeros(10,1);
%calculate first and second moment of all coefficients
%Initialize them here
b1mo=zeros(bigk,1);
b2mo=zeros(bigk,1);
%Number of burnin and kept draws
nburn=100;
nkeep=1000;
nrep=nburn+nkeep;
for irep=1:nrep
irep
%choose at random on of the bigk potential explanatory variables
%if it is already in the model, delete it else add it
%Based on this, make up candidate model
indch=round(bigk*rand);
xnew=xold;
mnewdraw=molddraw;
if indch>0
if molddraw(indch,1)==1
isum=0;
for i=1:indch
isum=isum+molddraw(i,1);
end
xnew = [xold(:,1:isum) xold(:,isum+2:kold)];
mnewdraw(indch,1)=0;
else
isum=0;
for i=1:indch
isum=isum+molddraw(i,1);
end
xnew = [xold(:,1:isum+1) xraw(:,indch) xold(:,isum+2:kold)];
mnewdraw(indch,1)=1;
end
end
knew=sum(mnewdraw)+1;
xtxinv=inv(xnew'*xnew);
ymy = y'*y - y'*xnew*xtxinv*xnew'*y;
lprobnew = .5*knew*log(g1) -.5*(n-1)*log(g2*ymy + g1*yty);
%Now decide whether to accept candidate draw
if log(rand) < (lprobnew - lprobold)
xold=xnew;
lprobold=lprobnew;
molddraw=mnewdraw;
kold=knew;
end
if irep>nburn
%If new drawn model better than current top 10, add it to list
for i=1:10
if lprobold>=lprobtop10(i,1)
if sum(abs(molddraw - top10mod(:,i)))<.09
break
end
if i<10
lprobtop10(i+1:10,1)=lprobtop10(i:9,1);
top10mod(:,i+1:10) = top10mod(:,i:9);
top10count(i+1:10,1)=top10count(i:9,1);
end
lprobtop10(i,1)=lprobold;
top10mod(:,i)=molddraw;
top10count(i,1)=0;
break
end
end
for i=1:10
temp1=sum(abs(molddraw-top10mod(:,i)));
if temp1<.01
top10count(i,1)=top10count(i,1)+1;
break
end
end
inccount = inccount + molddraw;
msize=msize + kold;
%calculating posterior properties of coefficients means
%we have to write out full posterior
Q1inv = (1+g0)*xold'*xold;
Q0inv=g0*xold'*xold;
Q1=inv(Q1inv);
b1= Q1*xold'*y;
vs2 = (y-xold*b1)'*(y-xold*b1) + b1'*Q0inv*b1;
bcov = (vs2/(n-2))*Q1;
%the next bit of this is awkward, needed to find out if variable is
%included in the model and, if so, to find out where it is
summer1=1;
for i=1:bigk
bc=zeros(1,kold);
if molddraw(i,1)==1
summer1=summer1+1;
bc(1,summer1)=1;
bmean=bc*b1;
bvar =bc*bcov*bc';
b1mo(i,1)=b1mo(i,1) + bmean;
b2mo(i,1)=b2mo(i,1) + (bvar+bmean^2);
end
end
end
end
inccount=inccount./nkeep;
'proportion of models visited containing each variable'
varcount=cumsum(ones(bigk,1));
[varcount inccount]
'mean number of regressors in models'
msize/nkeep
'Prob of top 10 models, analytical (after deleting rest of models)'
lt1=lprobtop10 - max(lprobtop10);
exp(lt1)/sum(exp(lt1))
'Prob of top 10 models, numerical (after deleting rest of models)'
top10count./sum(top10count)
'Prob of top 10 models out of total number of models, numerical'
sum(top10count)/nkeep
b1mo=b1mo./nkeep;
b2mo=b2mo./nkeep;
bsd=sqrt(b2mo-b1mo.^2);
'Posterior mean and stand dev of each coefficient'
[b1mo bsd]