# MC for stochastic volatility model # 03/06/2005 creastovol=function(n,beta2,sigma2,phi) { # Simulation of the data # t=1,....,n, y_t=beta*exp(alpha_t/2)*u_t # alpha_1 ~ N(0,sigma2/(1-phi^2)) # t=2....,n, alpha_t=phi*alphz_{t-1}+sigma*v_t # u_t ~ N(0,1) and v_t ~ N(0,1) alpha=rep(0,n) alpha[1]=rnorm(1,0,sqrt(sigma2/(1-phi^2))) for (i in 2:n) alpha[i]=rnorm(1,phi*alpha[i-1],sqrt(sigma2)) y=rnorm(n,0,sqrt(beta2*exp(alpha))) list(y=y,alpha=alpha) } #------------------------------------------------------------------------------# lposshep=function(y,alpha,phi,sigma2,beta2) { # Log-posterior distribution up to a multiplicative constant # f(y_1,....,y_n,alpha_1,....,alpha_n,beta2,phi,sigma2) n=length(y) 0.5*log(1-phi^2)-0.5*(n+1)*log(beta2)-0.5*(n+1)*log(sigma2)- 0.5/sigma2*((1-phi^2)*alpha[1]^2+sum((alpha[2:n]-phi*alpha[1:(n-1)])^2))-0.5*sum(alpha)- 0.5/beta2*sum(y^2*exp(-alpha)) } #------------------------------------------------------------------------------# updalphashep=function(y,alphap,phi,sigma2,beta2) { # Proposal distribution for updating the latent vector alpha n=length(y) alpha=alphap lopro=0 for (i in 1:n) { if (i==1) { mu=phi*alphap[2] sigp2=1/sigma2+0.5*y[1]^2/exp(-mu)/beta2 mean=mu/sigma2+0.5*y[1]^2*exp(-mu)*(1+mu)/beta2-0.5 mean=mean/sigp2 } if (i!=1 & i!=n) { mu=phi*(alpha[i-1]+alphap[i+1])/(1+phi^2) sigp2=(1+phi^2)/sigma2+0.5*y[i]^2*exp(-mu)/beta2 mean=mu*(1+phi^2)/sigma2+0.5*y[i]^2*exp(-mu)*(1+mu)/beta2-0.5 mean=mean/sigp2 } if (i==n) { mu=phi*alpha[n-1] sigp2=1/sigma2+0.5*y[n]^2/exp(-mu)/beta2 mean=mu/sigma2+0.5*y[n]^2*exp(-mu)*(1+mu)/beta2-0.5 mean=mean/sigp2 } alpha[i]=rnorm(1,mean,sqrt(1/sigp2)) lopro=lopro+dnorm(alpha[i],mean,sqrt(1/sigp2),log=T) } list(alpha=alpha,lopro=lopro) } #------------------------------------------------------------------------------# loproalphashep=function(alpha,y,alphap,phi,sigma2,beta2) { # Weight associated to the proposal distribution of the latent vector alpha n=length(y) sigp2=rep(0,n) mu=sigp2 mean=sigp2 mu[1]=phi*alphap[2] mu[2:(n-1)]=phi*(alpha[1:(n-2)]+alphap[3:n])/(1+phi^2) mu[n]=phi*alpha[n-1] sigp2[1]=1/sigma2+0.5*y[1]^2/exp(-mu[1])/beta2 sigp2[2:(n-1)]=(1+phi^2)/sigma2+0.5*y[2:(n-1)]^2*exp(-mu[2:(n-1)])/beta2 sigp2[n]=1/sigma2+0.5*y[n]^2/exp(-mu[n])/beta2 mean[1]=mu[1]/sigma2+0.5*y[1]^2*exp(-mu[1])*(1+mu[1])/beta2-0.5 mean[2:(n-1)]=mu[2:(n-1)]*(1+phi^2)/sigma2+0.5*y[2:(n-1)]^2*exp(-mu[2:(n-1)])*(1+mu[2:(n-1)])/beta2-0.5 mean[n]=mu[n]/sigma2+0.5*y[n]^2*exp(-mu[n])*(1+mu[n])/beta2-0.5 mean=mean/sigp2 sum(dnorm(alpha,mean,sqrt(1/sigp2),log=T)) } #------------------------------------------------------------------------------# updparshep=function(y,alpha,phip) { # Proposal (Gibbs type) distribution for updating the parameters: sigma2, phi and beta2 n=length(y) suffi1=0.5*sum(y^2*exp(-alpha)) beta2=1/rgamma(1,0.5*(n-1),rate=suffi1) lopro=-suffi1/beta2-(0.5*(n+1))*log(beta2)+0.5*(n-1)*log(suffi1) suffi2=0.5*((1-phip^2)*alpha[1]^2+sum((alpha[2:n]-phip*alpha[1:(n-1)])^2)) sigma2=1/rgamma(1,0.5*(n-1),rate=suffi2) lopro=lopro-suffi2/sigma2-(0.5*(n+1))*log(sigma2)+0.5*(n-1)*log(suffi2) suffi3=sum(alpha[2:(n-1)]^2) mean=sum(alpha[2:n]*alpha[1:(n-1)])/suffi3 var=sigma2/suffi3 u=runif(1) phi=mean+sqrt(var)*qnorm((u*(pnorm((1-mean)/sqrt(var))-pnorm((-1-mean)/sqrt(var)))+pnorm((-1-mean)/sqrt(var)))) lopro=lopro+dnorm(phi,mean,sqrt(var),log=T)-log(pnorm((1-mean)/sqrt(var))-pnorm((-1-mean)/sqrt(var))) list(beta2=beta2,sigma2=sigma2,phi=phi,lopro=lopro) } #------------------------------------------------------------------------------# loproparshep=function(beta2,sigma2,phi,y,alpha,phip) { # Weight associated to the proposal distribution of the parameters: sigma2, phi and beta2 n=length(y) suffi1=0.5*sum(y^2*exp(-alpha)) suffi2=0.5*((1-phip^2)*alpha[1]^2+sum((alpha[2:n]-phip*alpha[1:(n-1)])^2)) lopro=-suffi1/beta2-(0.5*(n+1))*log(beta2)+0.5*(n-1)*log(suffi1) lopro=lopro-suffi2/sigma2-(0.5*(n+1))*log(sigma2)+0.5*(n-1)*log(suffi2) suffi3=sum(alpha[2:(n-1)]^2) mean=sum(alpha[2:n]*alpha[1:(n-1)])/suffi3 var=sigma2/suffi3 lopro=lopro+dnorm(phi,mean,sqrt(var),log=T)-log(pnorm((1-mean)/sqrt(var))-pnorm((-1-mean)/sqrt(var))) lopro } #------------------------------------------------------------------------------# hmstovol=function(niter,y) { # Hasting-Metropolis algorithm # niter : number of iterations n=length(y) alpha=rep(0,n) phi=rep(0,niter) sigma2=rep(0,niter) beta2=rep(0,niter) phi[1]=2*runif(1)-1 sigma2[1]=abs(rnorm(1,0.5,1)) beta2[1]=abs(rnorm(1,0.5,1)) for (t in 1:(niter-1)) { for (i in 1:n) { #Gibbs parameters for alpha[i] if (i==1) { mu=phi[t]*alpha[2] sigp2=1/sigma2[t]+0.5*y[1]^2*exp(-mu)/beta2[t] mean=mu/sigma2[t]+0.5*y[1]^2*exp(-mu)*(1+mu)/beta2[t]-0.5 mean=mean/sigp2 } if (i!=1 & i!=n) { mu=phi[t]*(alpha[i-1]+alpha[i+1])/(1+phi[t]^2) sigp2=(1+phi[t]^2)/sigma2[t]+0.5*y[i]^2*exp(-mu)/beta2[t] mean=mu*(1+phi[t]^2)/sigma2[t]+0.5*y[i]^2*exp(-mu)*(1+mu)/beta2[t]-0.5 mean=mean/sigp2 } if (i==n) { mu=phi[t]*alpha[n-1] sigp2=1/sigma2[t]+0.5*y[n]^2*exp(-mu)/beta2[t] mean=mu/sigma2[t]+0.5*y[n]^2*exp(-mu)*(1+mu)/beta2[t]-0.5 mean=mean/sigp2 } prop=rnorm(1,mean,sqrt(1/sigp2)) if (i==1) { d1=exp(-0.5*(prop^2-2*phi[t]*prop*alpha[2])/sigma2[t]-0.5*prop-0.5*exp(-prop)*y[1]^2/beta2[t]) d2=exp(-0.5*(alpha[1]^2-2*phi[t]*alpha[2]*alpha[1])/sigma2[t]-0.5*alpha[1]-0.5*exp(-alpha[1])*y[1]^2/beta2[t]) } if (i!=1 & i!=n) { d1=exp(-0.5*(prop^2*(1+phi[t]^2)-2*phi[t]*prop*(alpha[i-1]+alpha[i+1]))/sigma2[t]-0.5*prop-0.5*exp(-prop)*y[i]^2/beta2[t]) d2=exp(-0.5*(alpha[i]^2*(1+phi[t]^2)-2*phi[t]*alpha[i]*(alpha[i-1]+alpha[i+1]))/sigma2[t]-0.5*alpha[i]-0.5*exp(-alpha[i])*y[i]^2/beta2[t]) } if (i==n) { d1=exp(-0.5*(prop^2-2*phi[t]*prop*alpha[n-1])/sigma2[t]-0.5*prop-0.5*exp(-prop)*y[n]^2/beta2[t]) d2=exp(-0.5*(alpha[n]^2-2*phi[t]*alpha[n-1]*alpha[n])/sigma2[t]-0.5*alpha[n]-0.5*exp(-alpha[n])*y[n]^2/beta2[t]) } prob=(d1*dnorm(alpha[i],mean,sqrt(1/sigp2)))/(d2*dnorm(prop,mean,sqrt(1/sigp2))) if (runif(1)