-
Notifications
You must be signed in to change notification settings - Fork 15
/
run_sp500_3parms_new_method.R
92 lines (76 loc) · 3.49 KB
/
run_sp500_3parms_new_method.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
#########################################################################
## Application on S&P500 using 3-parm SV Model ##
## 3 parameter model ##
## NEW joint sampling not adaptive ##
#########################################################################
source("R/pgasSV_binorm.R")
library(astsa)
library(mcmc)
library(plyr)
library(MASS)
load('data/sp500.gr.rda')
y = 100*window(sp500.gr, start=2005)
#-------------------------------------------------------------------
npart = 10 # Number of particles used in pgas
nmcmc = 2000 # Number of iterations in the mcmc samplers after burnin
burnin = 100 # Number of iterations to burn
mcmseed = 90210
# Bi-normal prior and initial values
init = c(.9, 0.3, .15)
hyper = c(0.9, 0.3, 0.1, 0.25, -0.6)
sigma_MH = .05*matrix(c(1,-0.25,-0.25,1),nrow=2,ncol=2)
# Run it
ptm <- proc.time()
u3na = pgasSV_binorm(nmcmc, burnin, y, init, hyper, sigma_MH, npart, parms_2 = FALSE, mcmseed)
(time2run = proc.time() - ptm)
# Acceptance Rate
cat("The acceptance rate is", u3na$acp, "%.")
### Pretty pictures
parms = cbind(u3na$phi, sqrt(u3na$q), u3na$mu)
names = c(expression(phi), expression(sigma), expression(mu))
culer = c(rgb(.66,.12,.85), rgb(.12,.66,.85), rgb(.8*.52,.8*.87,.8*.08) )
culerb = c(rgb(.66,.12,.85, .4), rgb(.12,.66,.85,.4), rgb(.8*.52,.8*.87,.8*.08,.4) )
# parameters
dev.new(width=9, height=6)
IF = sprintf("%.2f",
round(cbind(initseq(parms[,1])$var.pos/initseq(parms[,1])$gamma0,
initseq(parms[,2])$var.pos/initseq(parms[,2])$gamma0,
initseq(parms[,3])$var.pos/initseq(parms[,3])$gamma0), 3 ))
par(mfrow=c(3,3))
for (i in 1:3){
tsplot(parms[,i], ylab='trace', xlab='index', ,main='', col=culer[i])
abline(h=mean(parms[,i]), col=rgb(0,0,0,.5), lwd=1)
mtext(names[i], 3, line=.25, cex=1)
}
L=-.2; U=1
acf1(parms[,1], 100, ylim=c(L,U), main='', col=culer[1])
legend('topright', legend=paste('IF =', IF[1]), bty='n', cex=1.25)
acf1(parms[,2], 100, ylim=c(L,U), main='', col=culer[2])
legend('topright', legend=paste('IF =', IF[2]), bty='n', cex=1.25)
acf1(parms[,3], 100, ylim=c(L,U), main='', col=culer[3])
legend('topright', legend=paste('IF =', IF[3]), bty='n', cex=1.25)
for (i in 1:2){
hist(parms[,i], prob=TRUE, breaks=20, main='', xlab='')
Grid(nx=0, ny=NULL)
hist(parms[,i], breaks=20, col=culerb[i], border=culerb[i], prob=TRUE, add=TRUE)
lines(density(parms[,i], adjust=3))
abline(v=mean(parms[,i]), col=rgb(0,0,0,.8), lwd=2, lty=1)
}
hist(parms[,3], prob=TRUE, breaks=20, main='', xlab='', ylim=c(0,3.6), xlim=c(-1.5,1.5))
Grid(nx=0, ny=NULL)
hist(parms[,3], breaks=20, col=culerb[3], border=culerb[3], prob=TRUE, add=TRUE)
lines(density(parms[,3], adjust=4))
abline(v=mean(parms[,3]), col=rgb(0,0,0,.8), lwd=2, lty=1)
# states
tspar = tsp(sp500.gr)
mX = ts(apply(u3na$X, 2, mean), start=tspar[1], frequency=tspar[3])
lX = ts(apply(u3na$X, 2, quantile, 0.025), start=tspar[1], frequency=tspar[3])
uX = ts(apply(u3na$X, 2, quantile, 0.975), start=tspar[1], frequency=tspar[3])
culera = c(rgb(.85,.30,.12,.3), rgb(.12,.65,.85,.3))
dev.new(width=9, height=6)
par(mfrow=c(2,1))
tsplot(y, ylab='S&P500 Returns (%)', col='dodgerblue3', lwd=2)
tsplot(mX, col=culer[2], ylim=c(min(lX)-.1, max(uX)+.1), lwd=2, ylab='Log Volatiltiy (%)')
xx=c(time(mX), rev(time(mX)))
yy=c(lX, rev(uX))
polygon(xx, yy, border=NA, col=culera[2])