搜索
人大经济论坛 › 附件下载

附件下载

所在主题：
文件名:  joe_99m——hansen.zip
资料下载链接地址: https://bbs.pinggu.org/a-1018542.html
本附件包括：
  • data.m
  • data.out
  • DATA.PRG
  • graphic.tkf
  • invest.doc
  • INVEST.FMT
  • invest.txt
  • thresh.out
  • THRESH_P.asv
  • THRESH_P.m
  • THRESH_P.PRG
  • Readme.txt
  • invest.dta
附件大小:
737.76 KB   举报本内容
附件中是hansen的门槛回归的matlab程序和论文，
由于对matlab实在不了解，希望某位大神可以把附件中的matlab程序改写为R程序，
非常感谢！

如果能改写为类似lm(formula,data="")可以直接修改回归方程，读取变量和数据的形式，奖励可以增至200，
如果能修改至动态面板，奖励增加到500

非常感谢！！

这是matlab程序，方便没有装matlab的同学查看


  1. %%%%%%%%%%%%%%
  2. % THRESH_P.M %
  3. %%%%%%%%%%%%%%
  4. % This is a Matlab program file.
  5. % It replicates the estimation, testing and graphs reported in
  6. % "Threshold Effects in Non-Dynamic Panels:
  7. % Estimation, Testing and Inference"

  9. % For questions, please contact

  11. % Bruce E. Hansen
  12. % Department of Economics
  13. % Social Science Building
  14. % University of Wisconsin
  15. % Madison, WI 53706-1393
  16. % bhansen@ssc.wisc.edu
  17. % http://www.ssc.wisc.edu/~bhansen/


  20. % This program file loads the GAUSS dataset "invest.fmt".
  21. % It creates the output file "thresh.out"
  22. function thresh_p;
  23. global nt;
  24. global t;
  25. global n;
  26. global max_lag;
  27. global thresh;
  28. global cf;
  29. global qn;
  30. global qq1;
  31. global vgraph_;
  32. global tt;
  33. global yt;
  34. global xt;
  35. global ct;
  36. global ty;
  37. global cc;
  38. global k;
  39. global tt;
  40. global qn1;

  42. load invest.txt;
  43. t = 15;
  44. nt = length(invest(:,1));
  45. n = nt/t;

  47. i = invest(:,1); % investment/assets
  48. q = invest(:,2); % Tobin's Q
  49. c = invest(:,3); % cash-flow/assets
  50. d = invest(:,4); % debt/assets


  53. qn = 400; % number of quantiles to examine %
  54. conf_lev = .95; % confidence level for threshold %
  55. vgraph_ = 1; % set to 1 to graph likelihood ratios %
  56. boot_1_ = 300; % # of replications, 0 for no bootstrap, single (300) %
  57. boot_2_ = 300; % # of replications, 0 for no bootstrap, double (300) %
  58. boot_3_ = 300; % # of replications, 0 for no bootstrap, triple (300) %
  59. trim_1_ = .01; % percentage to trim before search, single %
  60. trim_2_ = .01; % percentage to trim before search, double %
  61. trim_3_ = .05; % percentage to trim before search, triple %

  63. max_lag = 1;
  64. tt = t-max_lag;
  65. ty = n*(t-max_lag-1);

  67. y= lag_v(i,0);yt = tr(y);
  68. cf = lag_v(c,1);ct = tr(cf);
  69. q1 = lag_v(q,1);
  70. d1 = lag_v(d,1); % set to threshold variable %

  72. x=[q1,q1.^2,q1.^3,d1,q1.*d1];
  73. k=length(x(1,:));
  74. xt=zeros(length(yt(:,1)),k);
  75. j=1;
  76. while j<=k
  77. xt(:,j)=tr(x(:,j));
  78. j=j+1;
  79. end;


  82. thresh=d1;
  83. dd=unique(thresh);
  84. qnt1=qn*trim_1_;
  85. sq=trim_1_:(1/qn):(trim_1_+(1/qn)*(qn-2*qnt1));
  86. qq1=dd(floor(sq*length(dd(:,1))));
  87. qn1=length(qq1(:,1));
  88. cc=-2*log(1-sqrt(conf_lev));
  89. out=fopen('out.txt','wt');
  90. fprintf(out,'Number of Firms: %u\n',n);
  91. fprintf(out,'Number of Years used: %u\n',tt);
  92. fprintf(out,'Total Observations: %u\n',ty);
  93. fprintf(out,'Number of Quantiles: %u\n',qn);
  94. fprintf(out,'Confidence Level: %f\n',conf_lev);
  95. fprintf(out,'\n');
  96. fprintf(out,'\n');
  97. fprintf(out,'*********************************************\n');
  98. fprintf(out,'\n');
  99. fprintf(out,'\n');
  100. sse0 = sse_calc(yt,[xt,ct]);
  101. fprintf(out,'Zero Threshold Model\n');
  102. fprintf(out,'Sum of Squared Errors: %f\n',sse0);
  103. fprintf(out,'\n');
  104. fprintf(out,'\n');
  105. fprintf(out,'*********************************************\n');
  106. fprintf(out,'\n');
  107. fprintf(out,'\n');

  109. fprintf(out,'Single Threshold Model\n');
  110. fprintf(out,'\n');
  111. fclose(out);
  112. rhat1=model(0,trim_1_,boot_1_,0);
  113. out=fopen('out.txt','a+');
  114. fprintf(out,'*********************************************\n');
  115. fprintf(out,'\n');
  116. fprintf(out,'\n');

  118. fprintf(out,'Double Threshold Model\n');
  119. fprintf(out,'Trimming Percentage: %f\n',trim_2_);
  120. fprintf(out,'\n');
  121. fprintf(out,'First Iteration\n');
  122. fclose(out);
  123. rhat2=model(rhat1,trim_2_,boot_2_,2);
  124. out=fopen('out.txt','a+');
  125. fprintf(out,'Second Iteration\n');
  126. fclose(out);
  127. rhat1=model(rhat2,trim_2_,0,1);
  128. out=fopen('out.txt','a+');
  129. fprintf(out,'\n');
  130. fprintf(out,'\n');
  131. fprintf(out,'*********************************************\n');
  132. fprintf(out,'\n');
  133. fprintf(out,'\n');

  135. fprintf(out,'Triple Threshold Model\n');
  136. fprintf(out,'Trimming Percentage: %f\n',trim_3_);
  137. fprintf(out,'\n');
  138. fclose(out);
  139. rhat3=model([rhat1;rhat2],trim_3_,boot_3_,3);
  140. out=fopen('out.txt','a+');
  141. fprintf(out,'\n');
  142. fprintf(out,'\n');
  143. fprintf(out,'*********************************************\n');
  144. fprintf(out,'\n');
  145. fprintf(out,'\n');
  146. fclose(out);
  147. % Functions %
  148. function r=tr(y);
  149. global n;
  150. global tt;
  151. yf=reshape(y',tt,n)';
  152. yfm=yf-mean(yf')'*ones(1,length(yf(1,:)));
  153. yfm=yfm(:,1:tt-1)';
  154. ind=0;
  155. for i=1:length(yfm(1,:))
  156. for j=1:length(yfm(:,1))
  157. if ind==0
  158. r=yfm(j,i);
  159. ind=1;
  160. else
  161. r=[r;yfm(j,i)];
  162. end;
  163. end;
  164. end;



  168. function ss=sse_calc(y,x)
  169. e=y-x*(y'/x')';
  170. ss=e'*e;

  172. function r=lag_v(x,lagn);
  173. global max_lag;
  174. global t;
  175. global n;
  176. y=reshape(x',t,n)';
  177. y=y(:,(1+max_lag-lagn):(t-lagn))';
  178. ind=0;
  179. for i=1:length(y(1,:))
  180. for j=1:length(y(:,1))
  181. if ind==0
  182. r=y(j,i);
  183. ind=1;
  184. else
  185. r=[r;y(j,i)];
  186. end;
  187. end;
  188. end;

  190. function sse=thr_sse(y,q,r);
  191. global thresh;
  192. global cf;
  193. global xt;
  194. global ct;
  195. n=length(q(:,1));
  196. sse=zeros(n,1);
  197. qi=1;
  198. while qi<=n
  199. if r==0
  200. rr=q(qi);
  201. else
  202. rr=[r;q(qi)];
  203. end;
  204. rr=sortrows(rr,1);
  205. xx=[xt,ct];
  206. j=1;
  207. while j<=length(rr(:,1));
  208. d=(thresh<rr(j));
  209. xx=[xx,tr(cf.*d)];
  210. j=j+1;
  211. end;
  212. sse(qi)=sse_calc(y,xx);
  213. qi=qi+1;
  214. end;

  216. function [rsse,rqq]=r_est(y,r,trim_);
  217. global qn;
  218. global qn1;
  219. global qq1;
  220. if max(r)'==0
  221. qq=qq1;
  222. rr=0;
  223. else
  224. rr=sortrows(r,1);
  225. i=(1:qn1)';
  226. nn=sum((qq1*ones(1,length(rr)))<(ones(length(qq1(:,1)),1)*(rr')));
  227. qnt=qn*trim_;
  228. ii=(((i*ones(1,length(nn)))<=(ones(length(i(:,1)),1)*(nn+qnt)))-((i*ones(1,length(nn)))<=(ones(length(i(:,1)),1)*(nn-qnt))))*ones(length(rr(:,1)),1);
  229. ind=0;
  230. for j=1:length(ii)
  231. if ii(j)==0
  232. if ind==0
  233. qq=qq1(j);
  234. ind=1;
  235. else
  236. qq=[qq;qq1(j)];
  237. end;
  238. end;
  239. end;
  240. end;
  241. sse=thr_sse(y,qq,rr);
  242. [temp,rihat]=min(sse);
  243. clear temp;
  244. rsse=sse(rihat);
  245. rqq=qq(rihat);

  247. function rhat=model(r,trim_,rep,it);
  248. global qq1;
  249. global vgraph_;
  250. global yt;
  251. global ty;
  252. global cc;
  253. global xt;
  254. global thresh;
  255. global cf;
  256. global n;
  257. global k;
  258. global ct;
  259. global tt;
  260. global qn1;
  261. global qn;
  262. if max(r)==0
  263. qq=qq1;
  264. rr=0;
  265. else
  266. rr=sortrows(r,1);
  267. i=(1:qn1)';
  268. ind=0;
  269. for ij=1:length(rr)
  270. if ind==0
  271. nn=sum(qq1<rr(ij));
  272. ind=1;
  273. else
  274. nn=[nn,sum(qq1<rr(ij))];
  275. end;
  276. end;
  277. qnt=qn*trim_;
  278. inn1=(i*ones(1,length(nn))<=(ones(length(i),1)*(nn+qnt)));
  279. inn2=(i*ones(1,length(nn))<=(ones(length(i),1)*(nn-qnt)));
  280. ii=(inn1-inn2)*ones(length(rr(:,1)),1);
  281. ind=0;
  282. for j=1:length(ii)
  283. if ii(j)==0
  284. if ind==0
  285. qq=qq1(j);
  286. ind=1;
  287. else
  288. qq=[qq;qq1(j)];
  289. end;
  290. end;
  291. end;
  292. end;
  293. sse=thr_sse(yt,qq,rr);
  294. [temp,rihat]=min(sse);
  295. clear temp;
  296. rhat=qq(rihat);
  297. sse1=sse(rihat);
  298. lr=(sse/sse1-1)*ty;
  299. ind=0;
  300. temp=(lr<cc);
  301. for j=1:length(qq(:,1))
  302. if temp(j)==1
  303. if ind==0
  304. rhats=qq(j);
  305. ind=1;
  306. else
  307. rhats=[rhats;qq(j)];
  308. end;
  309. end;
  310. end;
  311. if vgraph_==1
  312. figure;
  313. plot(qq,lr,qq,ones(length(qq),1)*cc);
  314. grid on;
  315. if it==0
  316. title('Figure 1 LConfidence Interval Construction in Single Threshold Model');
  317. xlabel('Threshold Parameter');
  318. end;
  319. if it==1
  320. title('Figure 3 LConfidence Interval Construction in Double Threshold Model');
  321. xlabel('First Threshold Parameter');
  322. end;
  323. if it==2
  324. title('Figure 2 LConfidence Interval Construction in Double Threshold Model');
  325. xlabel('Second Threshold Parameter');
  326. end;
  327. if it==3
  328. title('Confidence Interval Construction in Triple Threshold Model');
  329. xlabel('Third Threshold Parameter');
  330. end;
  331. ylabel('Likelihood Ratio');
  332. end;
  333. out=fopen('out.txt','a+');
  334. if abs(max(r)')>0
  335. fprintf(out,'Fixed Thresholds: \n');
  336. for i=1:length(rr)
  337. fprintf(out,'%f ',rr(i));
  338. end;
  339. fprintf(out,'\n');
  340. rrr=sortrows([rr;rhat],1);
  341. else
  342. rrr=rhat;
  343. end;

  345. fprintf(out,'Threshold Estimate: %f\n',rhat);
  346. fprintf(out,'Confidence Region: %f %f\n',min(rhats),max(rhats));
  347. fprintf(out,'Sum of Squared Errors: %f\n',sse1);
  348. fprintf(out,'Trimming Percentage: %f\n',trim_);
  349. fprintf(out,'\n');
  350. fprintf(out,'\n');

  352. nr=length(rrr(:,1));
  353. xx=xt;
  354. dd=zeros(length(thresh(:,1)),nr);
  355. j=1;
  356. while j<=nr
  357. dd(:,j)=(thresh<rrr(j));
  358. d=dd(:,j);
  359. if j>1;
  360. d=d-dd(:,j-1);
  361. end;
  362. xx=[xx,tr(cf.*d)];
  363. j=j+1;
  364. end;
  365. d=1-dd(:,nr);
  366. xx=[xx,tr(cf.*d)];
  367. xxi=inv(xx'*xx);
  368. beta=xxi*(xx'*yt);
  369. e=yt-xx*beta;
  370. xxe=xx.*(e*ones(1,length(xx(1,:))));
  371. sehet=sqrt(diag(xxi*xxe'*xxe*xxi));
  372. sehomo=sqrt(diag(xxi*(e'*e))/(ty-n-length(xx(1,:))));
  373. fprintf(out,'Thresholds');
  374. for j=1:length(rrr)
  375. fprintf(out,'%f ',rrr(j));
  376. end;
  377. fprintf(out,'\n');
  378. fprintf(out,'\n');
  379. fprintf(out,'Regime-independent Coefficients standard errorshet standard errors\n');
  380. for j=1:k
  381. fprintf(out,'%f %f %f\n',beta(j),sehomo(j),sehet(j));
  382. end;
  383. fprintf(out,'\n');
  384. fprintf(out,'Regime- dependent Coefficients standard errorshet standard errors\n');
  385. for j=k+1:k+nr+1
  386. fprintf(out,'%f %f %f\n',beta(j),sehomo(j),sehet(j));
  387. end;
  388. fprintf(out,'\n');
  389. fprintf(out,'\n');

  391. if rep>0
  392. xx=[xt,ct];
  393. if abs(max(rr))>0
  394. j=1;
  395. while j<=length(rr(:,1))
  396. xx=[xx,tr(cf.*(thresh<rr(j)))];
  397. j=j+1;
  398. end;
  399. end;
  400. yp=xx*(yt'/xx')';
  401. e=yt-yp;
  402. sse0=e'*e;
  403. lrt=(sse0/sse1-1)*ty;
  404. fprintf(out,'LR Test for threshold effect: %f\n',lrt);
  405. fprintf(out,'\n');
  406. fprintf(out,'\n');
  407. stats=zeros(rep,1);

  409. j=1;
  410. while j<=rep
  411. eb=reshape(e',tt-1,n)';
  412. ind=0;
  413. y=eb(ceil(unifrnd(0,1,n,1)*n),:)';
  414. for i=1:length(y(1,:))
  415. for jjj=1:length(y(:,1))
  416. if ind==0
  417. rrrrr=y(jjj,i);
  418. ind=1;
  419. else
  420. rrrrr=[rrrrr;y(jjj,i)];
  421. end;
  422. end;
  423. end;
  424. yb=yp+rrrrr;
  425. clear rrrrr;
  426. sse0=sse_calc(yb,[xt,ct]);
  427. [sse1,rhat_b]=r_est(yb,0,trim_);
  428. rrr=rhat_b;
  429. if abs(max(r))>0
  430. jj=1;
  431. while jj<=length(r(:,1))
  432. sse0=sse1;
  433. [sse1,rhat_b]=r_est(yb,rrr,trim_);
  434. rrr=[rrr;rhat_b];
  435. jj=jj+1;
  436. end;
  437. end;
  438. lrt_b=(sse0/sse1-1)*ty;
  439. stats(j)=lrt_b;
  440. fprintf(out,'Bootstrap Replication: %f %f\n',j,lrt_b);
  441. j=j+1;
  442. end;
  443. fprintf(out,'\n');
  444. fprintf(out,'\n');
  445. stats=sortrows(stats,1);
  446. crits=stats(ceil([.90;.95;.99]*rep));
  447. fprintf(out,'Number of Bootstrap replications:%f\n',rep);
  448. fprintf(out,'Bootstrap p-value : %f\n',mean(stats>lrt));
  449. fprintf(out,'Critical Values: %f\n',crits);
  450. fprintf(out,'\n');
  451. fprintf(out,'\n');
  452. fclose(out);
  453. end;
复制代码


    熟悉论坛请点击新手指南
下载说明
1、论坛支持迅雷和网际快车等p2p多线程软件下载，请在上面选择下载通道单击右健下载即可。
2、论坛会定期自动批量更新下载地址,所以请不要浪费时间盗链论坛资源,盗链地址会很快失效。
3、本站为非盈利性质的学术交流网站,鼓励和保护原创作品，拒绝未经版权人许可的上传行为。本站如接到版权人发出的合格侵权通知，将积极的采取必要措施；同时，本站也将在技术手段和能力范围内，履行版权保护的注意义务。
(如有侵权，欢迎举报)
二维码

扫码加我 拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

京ICP备16021002号-2 京B2-20170662号 京公网安备 11010802022788号 论坛法律顾问：王进律师 知识产权保护声明   免责及隐私声明

GMT+8, 2025-12-21 13:04