# STOCHASTIC VOLATILITY set.seed(1243) alpha = -0.00645 phi = 0.99 W = 0.15^2 n = 1000 h = rep(0,n) h[1] = rnorm(1,alpha/(1-phi),sqrt(W/(1-phi^2))) for (t in 2:n) h[t] = rnorm(1,alpha+phi*h[t-1],sqrt(W)) sig = exp(h/2) y = rnorm(n,0,sig) plot(y,type="l") # Bayesian sequential learning log-volatilities # via particle filter (bootstrap filter!) m0 = 0 C0 = 1 # Filtro de particulas set.seed(1234) M = 1000 h = rnorm(M,m0,sqrt(C0)) hs = matrix(0,M,n) for (t in 1:n){ # propagacao h1 = rnorm(M,alpha+phi*h,sqrt(W)) # resampling w = dnorm(y[t],0,exp(h1/2)) h = sample(h1,replace=TRUE,size=M,prob=w) hs[,t] = h } qstdev = t(apply(exp(hs/2),2,quantile,c(0.05,0.5,0.95))) par(mfrow=c(1,1)) ts.plot(qstdev,col=c(3,2,3)) lines(sig,col=4) par(mfrow=c(3,5)) for (t in trunc(seq(1,n,length=15))){ hist(exp(hs[,t]/2),main=paste("t=",t,sep=""),xlab="st.dev.",prob=TRUE) abline(v=sig[t],col=2) }