winbugs程序问题

epoh 发表于 2011-8-9 20:38
Stochastic volatility model:  1.Observation Equations  2.State Equations数学式子不好表达详细请参考sv ...

#SV-M AR(1) MODEL:
#y[t]=a+b*y[t-1]+d*exp(theta[t])+exp(theta[t]/2)*u[t],u[t]iid~N(0,1)
#theta[t]=mu+phi*(theta[t-1]-mu)+v[t], v[t]iid~N(0,1/itau2)

model
{
  beta1<-a/(1-b)
  beta2<-(1-b*b)/exp(theta0)
  y0~dnorm(beta1,beta2)
  ymean[1]<-a+b*y0+d*exp(theta[1])
  yisigma2[1]<-1/exp(theta[1])
  y[1]~dnorm(ymean[1],yisigma2[1])
  for(i in 2:n)
  {
   ymean<-a+b*y[i-1]+d*exp(theta)
   yisigma2<-1/exp(theta)
   y~dnorm(ymean,yisigma2)
  }
  
  theta0~dnorm(mu,isigma)
  thetamean[1]<-mu+phi*(theta0-mu)
  theta[1]~dnorm(thetamean[1],itau2)
  for(j in 2:n)
  {
    thetamean[j]<-mu+phi*(theta[j-1]-mu)
    theta[j]~dnorm(thetamean[j],itau2)
  }

  a~dnorm(1,1.0E-6)
  psi~dnorm(0,1.0E-6)  #the question is why here add psi variable
  b<-(exp(psi)-1)/(1+exp(psi))
  phi1~dbeta(20,1.5)
  phi<-2*phi1-1
  tau<-sqrt(1/itau2)
  mu~dnorm(0,0.01)
  itau2~dgamma(2.5,0.025)
  d~dnorm(0,1)
  isigma<-(1-phi*phi)*itau2
  
}

#06year rate
list(y=c(…………),n)

list(a=0,mu=0,phi1=0.95,itau2=50,d=0,psi=0)

以上可以通过，但我不知道代码是否对照了上面的模型，我的疑问有2：
1)为什么计算b时要加入psi变量？而不是直接b<-dnorm(1,1.0E-6)。如果是要加入的，那为什么计算公式是b<-(exp(psi)-1)/(1+exp(psi))
2)对于一般SV模型，当然这个也一样，theta0按理说应该是服从dnorm(mu,isigma),其中isigma<-(1-phi*phi)*itau2
  但是很多程序，包括你在一些帖子中提供的程序，都是theta0~dnorm(mu,itau2)
  我的发现是：
  根据theta[t]=mu+phi*(theta[t-1]-mu)+v[t]
  则有var(theta0)=phi*phi*var(theta0)+1/itau2
  推出var(theta0)=1/((1-phi*phi)*itau2),在程序中，倒一下。这是我从上面y0的方差计算中得到的计算方式，不知对否。
望老师指点下。谢谢。

openbugs是winbugs的开源版本 二者几乎一致 不会遇到关于version的we你

liuxin9023 发表于 2011-8-12 17:01
openbugs是winbugs的开源版本 二者几乎一致 不会遇到关于version的we你

你好，能帮我解决11楼的问题吗？呵呵

epoh 发表于 2011-8-9 20:38
Stochastic volatility model:  1.Observation Equations  2.State Equations数学式子不好表达详细请参考sv ...

#SV-M AR(1) MODEL:
#y[t]=a+b*y[t-1]+d*exp(theta[t])+exp(theta[t]/2)*u[t],u[t]iid~N(0,1)
#theta[t]=mu+phi*(theta[t-1]-mu)+v[t], v[t]iid~N(0,1/itau2)

model
{
  beta1<-a/(1-b)
  beta2<-(1-b*b)/exp(theta0)
  y0~dnorm(beta1,beta2)
  ymean[1]<-a+b*y0+d*exp(theta[1])
  yisigma2[1]<-1/exp(theta[1])
  y[1]~dnorm(ymean[1],yisigma2[1])
  for(i in 2:n)
  {
   ymean<-a+b*y[i-1]+d*exp(theta)
   yisigma2<-1/exp(theta)
   y~dnorm(ymean,yisigma2)
  }
  
  theta0~dnorm(mu,isigma)
  thetamean[1]<-mu+phi*(theta0-mu)
  theta[1]~dnorm(thetamean[1],itau2)
  for(j in 2:n)
  {
    thetamean[j]<-mu+phi*(theta[j-1]-mu)
    theta[j]~dnorm(thetamean[j],itau2)
  }

  a~dnorm(1,1.0E-6)
  psi~dnorm(0,1.0E-6)  #the question is why here add psi variable
  b<-(exp(psi)-1)/(1+exp(psi))
  phi1~dbeta(20,1.5)
  phi<-2*phi1-1
  tau<-sqrt(1/itau2)
  mu~dnorm(0,0.01)
  itau2~dgamma(2.5,0.025)
  d~dnorm(0,1)
  isigma<-(1-phi*phi)*itau2
  
}

#06year rate
list(y=c(…………),n)

list(a=0,mu=0,phi1=0.95,itau2=50,d=0,psi=0)

以上可以通过，但我不知道代码是否对照了上面的模型，我的疑问有2：
1)为什么计算b时要加入psi变量？而不是直接b<-dnorm(1,1.0E-6)。如果是要加入的，那为什么计算公式是b<-(exp(psi)-1)/(1+exp(psi))
2)对于一般SV模型，当然这个也一样，theta0按理说应该是服从dnorm(mu,isigma),其中isigma<-(1-phi*phi)*itau2
  但是很多程序，包括你在一些帖子中提供的程序，都是theta0~dnorm(mu,itau2)
  我的发现是：
  根据theta[t]=mu+phi*(theta[t-1]-mu)+v[t]
  则有var(theta0)=phi*phi*var(theta0)+1/itau2
  推出var(theta0)=1/((1-phi*phi)*itau2),在程序中，倒一下。这是我从上面y0的方差计算中得到的计算方式，不知对否。
望老师指点下。谢谢。
本文来自: 人大经济论坛 Matlab及其他计量软件专版 版，详细出处参考： https://bbs.pinggu.org/forum.php? ... amp;from^^uid=1569108

金币啊

ychaust 发表于 2011-10-18 21:39
金币啊

你回答啊。。。

有收获

#logistic regression for shorks
#breeding success as a function of both rainfall and temperature in May
model
{
  for(i in 1:N)
    {
     nbsuccess[i]~dbin(p[i],nbpairs[i])
     logit(p[i])<-intercept+slope.temperature*(temperature[i]-mean(temperature[]))/(sd(temperature[]))
       +slope.rainfall*(rainfall[i]-mean(rainfall[]))/(sd(rainfall[]))
    }
   #priors for regression parameters
   intercept~dnorm(0.0,0.001)
   slope.temperature~dnorm(0.0,0.001)
   slope.rainfall~dnorm(0.0,0.001)
}
#data
list(nbsuccess=c(151,105,73,107,113,87,77,108,118,122,112,120,122,89,69,71,53,
41,53,31,35,14,18),nbpairs=c(173,164,103,113,122,112,98,121,132,136,133,137,145,117,
90,80,67,54,58,39,42,23,23),temperature=c(15.1,13.3,15.3,13.3,14.6,15.6,13.1,13.1,15,11.7,15.3,14.4,14.4,12.7,11.7,11.9,15.9,13.4,14,13.9,12.9,15.1,13),
rainfall=c(67,52,88,61,32,36,72,43,92,32,86,28,57,55,66,26,28,96,48,90,86,78,87),
N=23)
#initial values
list(intercept=0,slope.temperature=0,slope.rainfall=0)

能否把这个PPT上传???