极大似然估计

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

哪位高手能帮我看看，帮我写个极大似然估计的程序，模型和似然函数见附件。因为是新手 ，所以不懂，希望高手能帮帮忙，万分感谢！

扫码加我 拉你入群

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

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

分享0 收藏0 回帖

关键词：极大似然估计 极大似然 似然估计 万分感谢 似然函数 程序 模型

你的附件下载还要收流量币。能贴出来吗？

TaskShare 发表于 2013-7-20 14:40
你的附件下载还要收流量币。能贴出来吗？

模型是关于短期利率随机波动的模型
具体形式是：

其中

=
是利率漂移的均值回复项， 是随机波动项，即方差。  服从标准正太分布。


我有利率的日数据，按照这个模型来做最大似然估计。




似然函数：  ,  , 跟上面的定义一样。

chenhao1101 发表于 2014-2-25 23:12
模型是关于短期利率随机波动的模型
具体形式是：

贴出来，看不到公式了，我设置的貌似不用论坛币就可以下载啊！

chenhao1101 发表于 2014-2-25 23:13
贴出来，看不到公式了，我设置的貌似不用论坛币就可以下载啊！

麻烦你了！谢谢

chenhao1101 发表于 2014-2-25 23:12
模型是关于短期利率随机波动的模型
具体形式是：

2013/7/15号的帖子，2014/2/25号还没有解决？

nuomin 发表于 2014-2-27 20:03
2013/7/15号的帖子，2014/2/25号还没有解决？

是的啊，一直没人解读！

chenhao1101 发表于 2014-3-2 22:48
是的啊，一直没人解读！

我试着算了一下，结果不是太好，不过总比没有强不是？

chenhao1101 发表于 2014-3-2 22:48
是的啊，一直没人解读！

你的模型设定里的vr部分不能保证是大0，为了计算能进行下去，我用vr=abs(...)代替了原来的vr。代码是用R语言写的，没有加各种的错误识别和规避的部分，所以，实际数据如果不事先剔除一些极端值的话会有错误发生。

1. ##generate time series
   
2. set.seed(83)
   
3. epsilon <- rnorm(1500,0,1)
   
4. rts <- epsilon
   
5. vrr <- numeric(1500)
   
6. for(i in 2:1500){
   
7.     ur <- 0.01+0.06*rts[i-1]+0.05*rts[i-1]^2-0.04*rts[i-1]^3+0.03*rts[i-1]^4-0.02*rts[i-1]^5+0.01/rts[i-1]
   
8.     vr <- sqrt(0.0005+0.0006*rts[i-1]+0.0007*rts[i-1]^2)*epsilon[i]
   
9.     rts[i]=ur+vr+rts[i-1]
   
10.     vrr[i] <- 0.0005+0.0006*rts[i-1]+0.0007*rts[i-1]^2
    
11. }
    
12. View(rts)
    
13. y <- rts[500:1500]
    
14. plot(y,type="l")
    
15. rm(i,ur,vr)
    
16. 
    
17. ##transform data
    
18. dy <- diff(y,lag=1,differences=1)
    
19. est_data <- as.data.frame(cbind(dy,rep(1,1000),y[-1],y[-1]^2,y[-1]^3,y[-1]^4,y[-1]^5,1/y[-1]))
    
20. dimnames(est_data)[[2]] <- c("deltay","slop","y1","y2","y3","y4","y5","y6")
    
21. View(est_data)
    
22. 
    
23. ##code mle function
    
24. 
    
25. attach(est_data)
    
26. mlefunc <- function(x){
    
27.     sapply(1:11,function(z) assign(paste("x",z,sep=""),x[z],pos=1))
    
28.     urf <- x11*slop+x1*y1+x2*y2+x3*y3+x4*y4+x5*y5+x6*y6
    
29.     vrf <- abs(x7*slop+x8*y1+x9*y1^x10)
    
30.     sum(-log(vrf)-((deltay-urf)^2)/ifelse(vrf<1e-8,1e-4,vrf))
    
31. }
    
32. 
    
33. mlegrfunc <- function(x){
    
34.     sapply(1:11,function(z) assign(paste("x",z,sep=""),x[z],pos=1))
    
35.     ugr <- x11*slop+x1*y1+x2*y2+x3*y3+x4*y4+x5*y5+x6*y6
    
36.     vgr <- abs(x7*slop+x8*y1+x9*y1^x10)
    
37.     c(sum(2*(deltay-ugr)*y1/vgr),
    
38.       sum(2*(deltay-ugr)*y2/vgr),
    
39.       sum(2*(deltay-ugr)*y3/vgr),
    
40.       sum(2*(deltay-ugr)*y4/vgr),
    
41.       sum(2*(deltay-ugr)*y5/vgr),
    
42.       sum(2*(deltay-ugr)*y6/vgr),
    
43.       sum(-1/ifelse(vgr<1e-8,1e-4,vgr)+(deltay-ugr)^2/ifelse(vgr^2<1e-8,1e-4,vgr^2)),
    
44.       sum(-y1/ifelse(vgr<1e-8,1e-4,vgr)+y1*(deltay-ugr)^2/ifelse(vgr^2<1e-8,1e-4,vgr^2)),
    
45.       sum(-y1^x10/ifelse(vgr<1e-8,1e-4,vgr)+y1^x10*(deltay-ugr)^2/ifelse(vgr^2<1e-8,1e-4,vgr^2)),
    
46.       sum(-x9*log(abs(y1))*y1^x10/ifelse(vgr<1e-8,1e-4,vgr)+-x9*log(abs(y1))*y1^x10*(deltay-ugr)^2/ifelse(vgr^2<1e-8,1e-4,vgr^2)),
    
47.       sum(2*(deltay-ugr)/vgr)
    
48.       )
    
49. 
    
50. }
    
51. 
    
52. ##optmize mlefunc
    
53. sann_out <- optim(rep(0.01,11),mlefunc,method="SANN")
    
54. sann_out
    
55. bfgs_out <- optim(rep(0.01,11),fn=mlefunc,gr=mlegrfunc,method="BFGS")
    
56. bfgs_out
    
57. lbfgsb_out <- optim(rep(0.01,11),fn=mlefunc,gr=mlegrfunc,method="L-BFGS-B")
    
58. lbfgsb_out
    
59. neldermead_out <- optim(rep(0.01,11),fn=mlefunc,gr=mlegrfunc,method="Nelder-Mead")
    
60. neldermead_out
    
61. cg_out <- optim(rep(0.01,11),fn=mlefunc,gr=mlegrfunc,method="CG")
    
62. cg_out
    
63. #brent_out <- optim(rep(0.01,11),fn=mlefunc,gr=mlegrfunc,method="Brent")
    
64. #brent_out
    
65. 
    
66. detach(est_data)
    
67. rm(list=ls(pattern="^x"))
    
68. save.image("archmle.RData")

复制代码

计算的结果是用Nelder-Mead方法较接近真实值。用梯度函数的optimize方法较其他方法估计结果稳定。
注意：估计结果和真实值不等，非常期待大神能给我指出不足之处

nuomin 发表于 2014-3-3 08:38
你的模型设定里的vr部分不能保证是大0，为了计算能进行下去，我用vr=abs(...)代替了原来的vr。代码是用R语 ...

谢谢啊，请问你这个是用r 软件做的么？还有计算的参数估计有没有什么t 统计量什么的！万分感谢啊