学习学习。。。。。

thanks

ding...........

看看

谢谢分享

看看哈谢谢

thanks！

看看是什么？

呵呵，真的还想学习R呢，楼主好人

yhao

reply

听说R很强大

学习学习。。。。。

网上研讨会,

好东西

############################################################################################
#
#  STOCHASTIC VOLATILITY MODEL
#
############################################################################################
#
# HEDIBERT FREITAS LOPES
# Associate Professor of Econometrics and Statistics
# The University of Chicago Booth School of Business
# 5807 South Woodlawn Avenue
# Chicago, Illinois, 60637
# Email : hlopes@ChicagoGSB.edu
# URL: http://faculty.chicagobooth.edu/hedibert.lopes/research/
#
############################################################################################
rm(list=ls())
source("sv-ar1-routines.R")

# Simulating the data
# -------------------
set.seed(1243)
mu   = -0.00645
phi  =  0.99
tau2 =  0.15^2
tau  = sqrt(tau2)
h0   = 0.0
n    = 500
h    = rep(0,n)
h[1] = rnorm(1,mu+phi*h0,tau)
for (t in 2:n)
  h[t] = rnorm(1,mu+phi*h[t-1],tau)
vol = exp(h)
sd  = sqrt(vol)
y   = rnorm(n,0,sd)

# True values
# -----------
ptr = c(mu,phi,tau2)
h0tr = h0
htr = h
vtr = vol
sdtr= sd

pdf(file="simulateddata.pdf",width=12,height=12)
par(mfrow=c(2,1))
plot(y,type="l",xlab="time",ylab="")
plot(vol,type="l",xlab="time",ylab="")
dev.off()

# Prior hyperparameters
# ---------------------
t
#################################################################################################
# Defining the mixture of 7 normals that approximate the log-chi-square with one degree of freedom
# ------------------------------------------------------------------------------------------------
set.seed(1576)
ms  = c(-11.40039,-5.24321,-9.83726,1.50746,-0.65098,0.52478,-2.35859)
s2s = c(  5.79596, 2.61369, 5.17950,0.16735, 0.64009,0.34023, 1.26261)
qs  = c(  0.00730, 0.10556, 0.00002,0.04395, 0.34001,0.24566, 0.25750)
mm  = sum(qs*ms)
vv  = sum(qs*(s2s+ms^2))-(mm)^2
MM  = 10000
x  = seq(-20,10,length=MM)
den = exp(-(x)/2)*exp(-0.5*exp(x))*exp(x)/sqrt(2*pi)
mix = rep(0,MM)
for (i in 1:MM)
  mix = sum(qs*dnorm(x,ms,sqrt(s2s)))
norm   = dnorm(x,mm,sqrt(vv))

pdf(file="logchi2.pdf",width=10,height=7)
par(mfrow=c(1,1))
plot(x,den,ylab="density",main="",type='n')
lines(x,den,col=1,lty=1,lwd=2)
lines(x,mix,col=2,lty=2,lwd=2)
legend(-15,0.2,legend=c("log chi^2_1","mixture of 7 normals"),lty=1:2,col=1:2,lwd=c(2,2),bty="n")
dev.off()

# Initial values and specific MCMC setup
# --------------------------------------
y1   = log(y^2)
h0   = h0tr
h    = htr
mu   = ptr[1]
phi  = ptr[2]
tau2 = ptr[3]
vh   = 0.1

set.seed(12578)
hs = NULL
ps = NULL
for (iter in 1:(M0+M)){
  h     = svu(y1,h,mu,phi,tau2,m0,C0)
  var   = 1/(iC0+phi^2/tau2)
  mean  = var*(iC0m0+phi*(h[1]-mu)/tau2)
  h0    = rnorm(1,mean,sqrt(var))
  X     = cbind(1,c(h0,h[1:(n-1)]))
  par   = fixedpar(h[1:n],X,theta0,iV0,nu0,s02)
  mu    = par[1]
  phi   = par[2]
  tau2  = par[3]
  ps    = rbind(ps,par)
  hs    = rbind(hs,h)
}
hs.mix = hs
ps.mix = ps
vol    = exp(hs[(M0+1):(M0+M),])
mvol   = apply(vol,2,median)
lvol   = apply(vol,2,quant05)
uvol   = apply(vol,2,quant95)

names = c("mu","phi","tau2")
pdf(file="par-mixture.pdf",width=10,height=7)
par(mfrow=c(3,3))
for (i in 1:3){
  plot(ps[(M0+1):(M0+M),i],type="l",xlab="iteration",ylab="",main=names)
  acf(ps[(M0+1):(M0+M),i],main="")
  hist(ps[(M0+1):(M0+M),i],xlab="",ylab="",main="")
}
dev.off()

pdf(file="vol-mixture.pdf",width=10,height=7)
par(mfrow=c(1,1))
plot(vtr,type="l",ylim=range(lvol,uvol,vtr),xlab="time",ylab="")
lines(lvol,col=2)
lines(uvol,col=2)
lines(mvol,col=4)
lines(vtr,col=1)
legend(250,7,legend=c("TRUE","POSTERIOR MEAN","90%CI"),col=c(1,4,2),lty=rep(1,3),lwd=rep(1,3),bty="n")
dev.off()

acfmix=NULL
for (t in 1:n)
  acfmix = rbind(acfmix,acf(hs.nor[(M0+1):(M0+M),t],plot=FALSE)$acf)
pdf(file="acf-mixture.pdf",width=10,height=7)
ts.plot(t(acfmix),xlab="lag",ylab="ACF",col=gray(0.8))
dev.off()

pdf(file="acf-rw-ind-normal-mixture.pdf",width=14,height=9)
L = min(acfrw,acfind,acfnor,acfmix)
U = max(acfrw,acfind,acfnor,acfmix)
par(mfrow=c(2,2))
ts.plot(t(acfrw),xlab="lag",ylab="ACF",ylim=c(L,U),main="RANDOM WALK",col=gray(0.8))
ts.plot(t(acfind),xlab="lag",ylab="ACF",ylim=c(L,U),main="INDEPENDENT",col=gray(0.8))
ts.plot(t(acfnor),xlab="lag",ylab="ACF",ylim=c(L,U),main="NORMAL",col=gray(0.8))
ts.plot(t(acfmix),xlab="lag",ylab="ACF",ylim=c(L,U),main="MIXTURE",col=gray(0.8))
dev.off()


# Comparing the four schemes
# --------------------------
mvol.rw   = apply(exp(hs.rw[(M0+1):(M0+M),]),2,mean)
mvol.ind = apply(exp(hs.ind[(M0+1):(M0+M),]),2,mean)
mvol.nor  = apply(exp(hs.nor[(M0+1):(M0+M),]),2,mean)
mvol.mix  = apply(exp(hs.mix[(M0+1):(M0+M),]),2,mean)
L = min(mvol.rw,mvol.ind,mvol.nor,mvol.mix,vtr)
U = max(mvol.rw,mvol.ind,mvol.nor,mvol.mix,vtr)
pdf(file="volatilities.pdf",width=14,height=9)
ts.plot(mvol.rw,xlab="time",ylab="",ylim=c(L,U),col=2)
lines(mvol.ind,col=3)
lines(mvol.nor,col=4)
lines(mvol.mix,col=5)
lines(vtr,col=1)
legend(300,4,legend=c("TRUE","RANDOM WALK","INDEPENDENT","NORMAL","MIXTURE"),lty=rep(1,5),bty="n",col=1:5)
dev.off()

pdf(file="parameters.pdf",width=12,height=6)
par(mfrow=c(1,3))
plot(density(ps.rw[,1]),xlab="",ylab="Density",main=expression(mu),col=2,ylim=c(0,70))
lines(density(ps.ind[,1]),col=3)
lines(density(ps.nor[,1]),col=4)
lines(density(ps.mix[,1]),col=5)
abline(v=ptr[1],lwd=2)

plot(density(ps.rw[,2]),xlab="",ylab="Density",main=expression(phi),col=2,ylim=c(0,70))
lines(density(ps.ind[,2]),col=3)
lines(density(ps.nor[,2]),col=4)
lines(density(ps.mix[,2]),col=5)
abline(v=ptr[2],lwd=2)

plot(density(ps.rw[,3]),xlab="",ylab="Density",main=expression(tau^2),col=2,ylim=c(0,70))
lines(density(ps.ind[,3]),col=3)
lines(density(ps.nor[,3]),col=4)
lines(density(ps.mix[,3]),col=5)
abline(v=ptr[3],lwd=2)
dev.off()