最高悬赏500论坛币，帮忙改个matlab程序

附件中是hansen的门槛回归的matlab程序和论文，
由于对matlab实在不了解，希望某位大神可以把附件中的matlab程序改写为R程序，
非常感谢！

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

非常感谢！！

这是matlab程序，方便没有装matlab的同学查看
paper.pdf (286.74 KB)

2011-12-22 09:40:37 上传

joe_99m——hansen.zip (737.76 KB) 本附件包括：

• 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

2011-12-22 09:38:34 上传

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"
   
8. 
   
9. % For questions, please contact
   
10. 
    
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/
    
18. 
    
19. 
    
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;
    
41. 
    
42. load invest.txt;
    
43. t = 15;
    
44. nt = length(invest(:,1));
    
45. n = nt/t;
    
46. 
    
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              
    
51. 
    
52. 
    
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      %
    
62. 
    
63. max_lag = 1;
    
64. tt = t-max_lag;
    
65. ty = n*(t-max_lag-1);
    
66. 
    
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 %
    
71. 
    
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;
    
80. 
    
81. 
    
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');
     
108. 
     
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');
     
117. 
     
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');
     
134. 
     
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;
     
165. 
     
166. 
     
167. 
     
168. function ss=sse_calc(y,x)
     
169. e=y-x*(y'/x')';
     
170. ss=e'*e;
     
171. 
     
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;
     
189. 
     
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;
     
215. 
     
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);
     
246. 
     
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;
     
344. 
     
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');
     
351. 
     
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 errors  het 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 errors  het 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');
     
390. 
     
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);
     
408. 
     
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;

复制代码

老兄不用花钱的
这原本就有三种code
Gauss,matlab,R
http://www.ssc.wisc.edu/~bhansen/progs/joe_99.html
joe_99r
    joe_99r.zip (124.09 KB) 本附件包括：

• invest.dat
• thresh_p.R
• invest.doc
• DATA.R

2011-12-22 10:07:51 上传

epoh 发表于 2011-12-22 10:06
老兄不用花钱的
这原本就有三种code
Gauss,matlab,R

呵呵，epoh兄每次回复都那么及时，让我等真是不知如何感谢

另外，以前去hansen网站浏览过，没见过hansen放过R程序呀，呵呵，可能是自己太孤陋寡闻了~

看来R的确有成为主流计量软件的趋势了，呵呵

最后，感谢epoh兄了哈！！

奇怪，为什么我不能设置最佳答案呢，时间不到？

呵呵，看了下更新时间，关于R的程序都是2011年后更新的，看来hansen老师也对R有研究啦

补充一下，如果有其他同学能改成类似lm(formula,data="")等可以直接使用的函数，奖励依旧不变

可以设置出售帖，我来购买

呵呵

kk22boy 发表于 2011-12-22 10:41
呵呵，看了下更新时间，关于R的程序都是2011年后更新的，看来hansen老师也对R有研究啦

补充一下，如果有 ...

补个winrats code
joe_99_rats
   joe_99_rats.rar (120.56 KB)

2011-12-22 10:48:18 上传

Panel Regression - Estimation by Fixed Effects
Dependent Variable INVA
Panel(15) of Annual Data From   1//1974:01 To    565//1987:01
Usable Observations   7910      Degrees of Freedom  7339
Total Observations   8474      Skipped/Missing      564
Centered R**2     0.410889      R Bar **2   0.365135
Uncentered R**2   0.807075      T x R**2    6383.960
Mean of Dependent Variable      0.0887200480
Std Error of Dependent Variable 0.0619148484
Standard Error of Estimate      0.0493327859
Sum of Squared Residuals        17.861098726
Regression F(570,7339)                8.9803
Significance Level of F           0.00000000
Log Likelihood                   12875.03147

   Variable                  Coeff     Std Error       T-Stat        Signif
******************************************************************
1.  Q{1}                 0.0104   8.9317e-04     11.63377  0.00000000
2.  Q2{1}         2.1291e-04   2.5598e-05     -8.31758  0.00000000
3.  Q3{1}         1.1672e-06   1.9501e-07      5.98531  0.00000000
4.  DEBTA{1}         -0.0222   4.2453e-03     -5.23262  0.00000017
5.  QD{1}         1.6349e-03   1.4222e-03      1.14954  0.25036936
6.  CFA{1}               0.0715   4.5259e-03     15.79410  0.00000000

Single Threshold       0.01582
Sum of Squared Residuals      17.78551
F Test vs 0 Threshold      33.61933
Double Threshold       0.01582       0.54180
Sum of Squared Residuals      17.72983
F Test vs 1 Threshold      24.83754

epoh 发表于 2011-12-22 09:46
老兄不用花钱的
这原本就有三种code
Gauss,matlab,R

瞧这人品，刚刚地

epoh 发表于 2011-12-22 09:46
老兄不用花钱的
这原本就有三种code
Gauss,matlab,R

不知道这位仁兄对这个程序熟么 我最近要用这个程序  但得改函数形式  身为门外汉不知从哪下手 不知能否指点一二？ 非常感谢

好赞的东西，可惜才看到，谢谢楼主

非常好的东西，留个印迹。