22# Example of particle Metropolis-Hastings in a stochastic volatility model
33# The effect on mixing while varying N.
44#
5- # Johan Dahlin <liu (at) johandahlin.com.nospam>
5+ # Johan Dahlin <uni (at) johandahlin.com.nospam>
66# Documentation at https://github.com/compops/pmh-tutorial
77# Published under GNU General Public License
88# #############################################################################
@@ -58,19 +58,19 @@ if (!loadSavedWorkspace) {
5858 for (k in 1 : length(noParticles )) {
5959 # Save the current time
6060 ptm <- proc.time()
61-
61+
6262 for (i in 1 : noSimulations ) {
6363 # Run the particle filter
6464 res <- particleFilterSVmodel(y , theta , noParticles [k ])
65-
65+
6666 # Save the log-Likelihood estimate
6767 logLikelihoodEstimates [k , i ] <- res $ logLikelihood
6868 }
69-
69+
7070 # Compute the variance of the log-likelihood and computational time per sample
7171 logLikelihoodVariance [k ] <- var(logLikelihoodEstimates [k , ])
7272 computationalTimePerSample [k ] <- (proc.time() - ptm )[3 ] / noSimulations
73-
73+
7474 # Print to screen
7575 print(paste(paste(paste(paste(" Simulation: " , k , sep = " " ), " of " , sep = " " ), length(noParticles ), sep = " " ), " completed." , sep = " " ))
7676 print(paste(paste(paste(paste(" No. particles: " , noParticles [k ], sep = " " ), " requires " , sep = " " ), computationalTimePerSample [k ], sep = " " ), " seconds for computing one sample." , sep = " " ))
@@ -104,21 +104,21 @@ if (loadSavedWorkspace) {
104104 resTheta <- array (0 , dim = c(length(noParticles ), noIterations - noBurnInIterations + 1 , 3 ))
105105 computationalTimePerIteration <- rep(0 , length(noParticles ))
106106 acceptProbability <- rep(0 , length(noParticles ))
107-
107+
108108 for (k in 1 : length(noParticles )) {
109109 # Save the current time
110110 ptm <- proc.time()
111-
111+
112112 # Run the PMH algorithm
113113 res <- particleMetropolisHastingsSVmodel(y , initialTheta , noParticles [k ], noIterations , stepSize = proposals [k , ,])
114-
114+
115115 # Save the parameter trace
116116 resTheta [k , ,] <- res $ theta [noBurnInIterations : noIterations ,]
117-
117+
118118 # Compute acceptance probability and computational time per sample
119- computationalTimePerIteration [k ] <- (proc.time() - ptm )[3 ] / noIterations
119+ computationalTimePerIteration [k ] <- (proc.time() - ptm )[3 ] / noIterations
120120 acceptProbability [k ] <- mean(res $ proposedThetaAccepted [noBurnInIterations : noIterations ])
121-
121+
122122 # Print to screen
123123 print(paste(paste(paste(paste(" Simulation: " , k , sep = " " ), " of " , sep = " " ), length(noParticles ), sep = " " ), " completed." , sep = " " ))
124124 }
@@ -134,7 +134,7 @@ for (k in 1:length(noParticles)) {
134134 acf_mu <- acf(resTheta [k , , 1 ], plot = FALSE , lag.max = 250 )
135135 acf_phi <- acf(resTheta [k , , 2 ], plot = FALSE , lag.max = 250 )
136136 acf_sigmav <- acf(resTheta [k , , 3 ], plot = FALSE , lag.max = 250 )
137-
137+
138138 resThetaIACT [k , ] <- 1 + 2 * c(sum(acf_mu $ acf ), sum(acf_phi $ acf ), sum(acf_sigmav $ acf ))
139139 resThetaIACTperSecond [k , ] <- resThetaIACT [k , ] / computationalTimePerIteration [k ]
140140}
@@ -148,7 +148,7 @@ print(table)
148148# #############################################################################
149149if (tuneProposals ) {
150150 proposals <- array (0 , dim = c(length(noParticles ), 3 , 3 ))
151-
151+
152152 for (k in 1 : length(noParticles )) {
153153 proposals [k , , ] <- cov(resTheta [k , , ]) * 2.562 ^ 2 / 3
154154 }
0 commit comments