学科带头人
可计算经济学与计量经济学实证十年爱好者
- 阅读权限
- 255
- 威望
- 0 级
- 论坛币
- 386113 个
- 通用积分
- 413.4972
- 学术水平
- 127 点
- 热心指数
- 140 点
- 信用等级
- 103 点
- 经验
- 46944 点
- 帖子
- 1767
- 精华
- 0
- 在线时间
- 2466 小时
- 注册时间
- 2007-11-5
- 最后登录
- 2024-4-15
|
代码文件,他的是一个章节,不太清楚完整不;
复制下面的代码,保存为matlab 的脚本文档,暂且命名为 Demo_BMS - %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]
复制代码
|
-
总评分: 学术水平 + 5
热心指数 + 5
信用等级 + 5
查看全部评分
|