epoh，帮我看下结果，好像说不通啊

宝笙宝姑娘 发表于 2011-12-25 16:00
Rmat=matrix(0,4,13)
Rmat[1,7]=Rmat[2,8]=Rmat[3,11]=Rmat[4,12]=1
rvec=rep(0,4)

module("finmetrics")
hp.ibm = seriesMerge(hp.s, ibm.s)
hp.ibm=100*hp.ibm
hp.ibm.bekk = mgarch(hp.ibm~1, ~bekk(1,1))
hp.ibm.bekk

#####
Rmat=matrix(0,4,13)
Rmat[1,7]=Rmat[2,8]=Rmat[3,11]=Rmat[4,12]=1
Rmat
rvec=rep(0,4)
bHat=as.vector(coef(hp.ibm.bekk ))
avar=Rmat%*%vcov(hp.ibm.bekk )%*%t(Rmat)
Wald=t((Rmat%*%bHat-rvec))%*%solve(avar)%*%(Rmat%*%bHat-rvec)
Wald   #39.53947
#A 2 B
Rmat=matrix(0,2,13)
Rmat[1,7]=Rmat[2,11]=1
Rmat
rvec=rep(0,2)
bHat=as.vector(coef(hp.ibm.bekk ))
avar=Rmat%*%vcov(hp.ibm.bekk )%*%t(Rmat)
Wald=t((Rmat%*%bHat-rvec))%*%solve(avar)%*%(Rmat%*%bHat-rvec)
Wald    #18.1506

#B 2 A

Rmat=matrix(0,2,13)
Rmat[1,8]=Rmat[2,12]=1
Rmat
rvec=rep(0,2)
bHat=as.vector(coef(hp.ibm.bekk ))
avar=Rmat%*%vcov(hp.ibm.bekk )%*%t(Rmat)
Wald=t((Rmat%*%bHat-rvec))%*%solve(avar)%*%(Rmat%*%bHat-rvec)
Wald   #15.10734

> Wald
         [,1]
[1,] 8.948605

epoh,是我自己前面弄错了。这个数值不大。
因为我前面每次会按照你给的程序输入
as.numeric(wald.stat)
p.wald = 1 - pchisq(wald.stat,4)

然后就会出来一个相同的数字。
as.numeric(wald.stat)
p.wald = 1 - pchisq(wald.stat,4)
这两句代码到底是什么意思？

[,1]
[1,] 8.948605
>

epoh,这到底是为什么呀？为什么我的结果两个这么大，一个这么小？
难道我做的数据一点意义也没有吗？

宝笙宝姑娘 发表于 2011-12-25 16:20
Rmat=matrix(0,4,13)
> Rmat[1,7]=Rmat[2,8]=Rmat[3,11]=Rmat[4,12]=1
> rvec=rep(0,4)

的确,你的数据在不同软件
都得到相同结果.
要不然你的数据要先处理成
类似matlab的要求试试
data     - A t by k matrix of zero mean residuals

epoh 发表于 2011-12-25 18:01
的确,你的数据在不同软件
都得到相同结果.
要不然你的数据要先处理成

好的，我再想办法。
那么LR和Wald的P值到底怎么得到呢？
我觉得下面这个命令不对

as.numeric(wald.stat)
p.wald = 1 - pchisq(wald.stat,4)
因为只要用这个命令，就一定得到0

宝笙宝姑娘 发表于 2011-12-25 18:08
好的，我再想办法。
那么LR和Wald的P值到底怎么得到呢？
我觉得下面这个命令不对

as.numeric(wald.stat)
p.wald = 1 - pchisq(wald.stat,4)
这两句是对的,wald.stat在正常范围都会有值
限制一个系数
  p.wald =1 - pchisq(wald.stat,1)
限制二个系数
  p.wald =1 - pchisq(wald.stat,2)
限制四个系数
  p.wald =1 - pchisq(wald.stat,4)
#######
数据可能出在rhl,请你仔细检查
看图形就知道.
hp.ibm
  
rsz.rhl
  

epoh 发表于 2011-12-25 19:53
as.numeric(wald.stat)
p.wald = 1 - pchisq(wald.stat,4)
这两句是对的,wald.stat在正常范围都会有值 ...

我知道了。验证过以后确实可以得到P值。
非常感谢你一直为我解答疑问！

能不能请教下在winrat里面如何计算lr检验呢，

请教各位：
关于极值理论POT模型参数估计，即估计广义帕累托分布（GPD）参数，方法采用贝叶斯估计。Openbugs(Winbugs)程序如下：
   model
        {
                for( i in 1 : N )
                   {
                                x[i] ~ dgpar(mu, tau,eta)
                   }
               
        # Prior distributions of the model parameters       

        # Gamma Prior               
                        eta~ dpar(0.001, 0.001)
                                tau~ dgamma(0.001, 0.001)       
                                       
        }
        list(N=30,mu=62, x=c(63.3,67,75,77,80,90,100,100,100,100,100,103,120,136,140,145,150,150,190,270,300,300,400,400,500,700,800,827,1825,2400))
list(tau=1,eta=0.4)
运行结果
                mean        sd        MC_error        val2.5pc        median        val97.5pc        start        sample
        eta        0.897        0.3875        0.002003        0.2862        0.8464        1.796        20000        80001
        tau        106.9        41.31        0.2351        46.51        100.3        206.5        20000        80001
可见：参数tau标准差为41.31，很大。
另外，我对数据用matlab拟合广义帕累托分布，得到参数估计（最大似然法）的标准差也很大。我本意是希望贝叶斯估计比最大似然估计效果要好些，结果均不理想。
问题：
1.程序是否正确？利用该方法，结果如何改进？
2.想用Bootstrap方法结合贝叶斯估计，或单独用Bootstrap方法进行参数估计，我没弄过。
3.其他能改进结果的简单方法。
   

你解决了吗  能不能帮帮我啊