-
Notifications
You must be signed in to change notification settings - Fork 16
/
run_sp500_2parms_new_method.R
133 lines (105 loc) · 4.38 KB
/
run_sp500_2parms_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
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
##############################################
## SP500 SV Model using PGAS_adaptive ##
## 2 parameters model ##
## this is using the new method ##
##############################################
library(astsa)
source("R/pgasSV_binorm.R")
library(plyr)
library(MASS)
library(mcmc)
load('data/sp500.gr.rda')
y = 100*window(sp500.gr, start=2005)
#-------------------------------------------------------------------
npart = 20 # Number of particles used in pgas
nmcmc = 5000 # 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.35, 0)
hyper = c(0.9, 0.35, 0.125, 0.25, -0.6)
sigma_MH = .1*matrix(c(1,-0.25,-0.25,1),nrow=2,ncol=2)
# Run it
ptm <- proc.time()
uu = pgasSV_binorm(nmcmc, burnin, y, init, hyper, sigma_MH, npart, parms_2 = T, mcmseed)
(time2run = proc.time() - ptm)
# Acceptance Rate
cat("The acceptance rate is", uu$acp, "%.")
############### Results ########################
culerb = c(rgb(.85,.30,.12,.3), rgb(.12,.65,.85,.3))
culera = c(rgb(.85,.30,.12,.6), rgb(.12,.65,.85,.7))
culer = c(rgb(.85,.30,.12), rgb(.12,.65,.85))
### Pretty pictures
parms = cbind(uu$phi, sqrt(uu$q))
names = c(expression(phi), expression(sigma))
nobs = length(y)
tspar = tsp(sp500.gr)
mX = ts(apply(uu$X, 2, mean), start=tspar[1], frequency=tspar[3])
lX = ts(apply(uu$X, 2, quantile, 0.025), start=tspar[1], frequency=tspar[3])
uX = ts(apply(uu$X, 2, quantile, 0.975), start=tspar[1], frequency=tspar[3])
# state and data
dev.new(width=9, height=6)
par(mfrow=c(2,1))
tsplot(y, ylab='S&P500 Returns (%)', col='dodgerblue3', lwd=2)
tsplot(mX, col=culera[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=culerb[2])
# 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), 2))
layout(matrix(c(1,3, 2,4), 2), heights=c(1,1))
for (i in 1:2){
tsplot(parms[,i], ylab='trace', xlab='index', ,main='', col=culer[i])
abline(h=mean(parms[,i]), col=culera[i%%2+1], lwd=2)
mtext(names[i], 3, line=.25, cex=1)
}
u1 = acf1(parms[,1], 100, plot=FALSE)
u2 = acf1(parms[,2], 100, plot=FALSE)
L = min(u1,u2)
U = max(c(u1,u2,1))
tsplot(u2, col=culera[2], lwd=2, ylab='ACF', xlab='LAG', ylim=c(L,U))
lines(u1, col=culera[1], lwd=2, lty=6)
abline(h=0)
legend('topright', legend=c(as.expression(bquote(~~.(IF[1])~~"["~phi~"]")),
as.expression(bquote(.(IF[2])~~"["~sigma~"]"))),
lwd=2, col=culer, title.col=1, lty=c(6,1), text.col=culer, title="Inefficiency", bg='white', text.font=2)
plot(parms[,1], parms[,2], pch=20, col=rgb(24,116,205,max=255,alpha=150), xlab=expression(phi),
ylab=expression(sigma), cex=1.25, cex.lab=1.25, panel.first=Grid() )
abline(h=mean(parms[,2]), col=gray(.5))
abline(v=mean(parms[,1]), col=gray(.5))
# > cor(parms)
# [,1] [,2]
# [1,] 1.0000000 -0.5451042
# [2,] -0.5451042 1.0000000
# > colMeans(parms)
# [1] 0.8483645 0.4533468
############### data and states from new and old
tspar = tsp(sp500.gr)
# old stuff
mXo = ts(apply(u$X, 2, mean), start=tspar[1], frequency=tspar[3])
lXo = ts(apply(u$X, 2, quantile, 0.025), start=tspar[1], frequency=tspar[3])
uXo = ts(apply(u$X, 2, quantile, 0.975), start=tspar[1], frequency=tspar[3])
#####
# pdf(file="sp500_states_data.pdf", width=9, height=6)
dev.new(width=9, height=6)
culer = c(rgb(.85,.30,.12), rgb(.12,.65,.85))
culerb = c(rgb(.85,.30,.12,.3), rgb(.12,.65,.85,.3))
par(mfrow=c(3,1))
tsplot(y, ylab='S&P500 Returns (%)', col='dodgerblue3', lwd=2)
tsplot(mXo, col=culer[1], ylim=c(min(lXo)-.1, max(uXo)+.1), lwd=2, ylab='Log Volatiltiy', )
xx=c(time(mXo), rev(time(mXo)))
yy=c(lX, rev(uXo))
polygon(xx, yy, border=NA, col=culerb[1])
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=culerb[2])
################## remap parms to causal #################
# parme = parms
# nn = nrow(parms)
# for (i in 1:nn){
# if (parms[i,1]>1) {parme[i,2]= parms[i,2]/parms[i,1]; parme[i,1]=1/parms[i,1] }
# }