Updated code and uploaded saved workspace for R

Former-commit-id: 6f549f7 [formerly 6f549f7 [formerly 15c47c9]]
Former-commit-id: 08f8258094ac504c4a8d9b95fee14d991a460c7f
Former-commit-id: d3cd8bf

Author: Johan Dahlin

README.md

@@ -18,9 +18,11 @@ Included files
 
 **example2-lgss.[R,py,m]** Implements the numerical illustration in Section 3.2 of parameter estimation of the parameter phi in the LGSS model using particle Metropolis-Hastings (PMH) with the faPF as the likelihood estimator. The output is the estimated parameter posterior as presented in Figure 4.
 
-**example3-sv.[R,py,m]** Implements the numerical illustration in Section 4.1 of parameter estimation of the three parameters in the stochastic volatility (SV) model using particle Metropolis-Hastings (PMH) with the bootstrap particle filter (bPF) as the likelihood estimator. The output is the estimated parameter posterior as presented in Figure 5.
+**example3-sv.[R,py,m]** Implements the numerical illustration in Section 4.1 of parameter estimation of the three parameters in the stochastic volatility (SV) model using particle Metropolis-Hastings (PMH) with the bootstrap particle filter (bPF) as the likelihood estimator. The output is the estimated parameter posterior as presented in Figure 5. The code takes some time (hours to execute).
 
-**example3-sv.[R,py,m]** Implements the numerical illustration in Section 4.2, which makes use of the same setup as in Section 4.1 but with a tuned proposal distribution. The output is the estimated ACF and IACT as presented in Figure 7.
+**example4-sv.[R,py,m]** Implements the numerical illustration in Section 4.2, which makes use of the same setup as in Section 4.1 but with a tuned proposal distribution. The output is the estimated ACF and IACT as presented in Figure 7. The code takes some time (hours to execute).
+
+**example*.[RData,mat]** A saved copy of the workspace after running to code.
 
 Supporting files
 --------------

python/example1-lgss.py

@@ -0,0 +1,78 @@
+##############################################################################
+#
+# Example of particle filtering
+# in a linear Gaussian state space model
+#
+# (c) 2015 Johan Dahlin
+# johan.dahlin (at) liu.se
+#
+##############################################################################
+
+# Import some libraries
+import matplotlib.pylab as plt
+import numpy as np;
+import stateEstimationHelper as helpState
+
+# Set the random seed to replicate results in tutorial
+np.random.seed( 10 );
+
+##############################################################################
+# Define the model
+##############################################################################
+
+# Here, we use the following model
+#
+# x[tt+1] = phi   * x[tt] + sigmav * v[tt]
+# y[tt]   = x[tt]         + sigmae * e[tt]
+#
+# where v[tt] ~ N(0,1) and e[tt] ~ N(0,1)
+
+# Set the parameters of the model (par[0],par[1],par[2])
+par = np.zeros(3)
+par[0] = 0.75
+par[1] = 1.00
+par[2] = 1.00
+
+# Set the number of time steps to simulate
+T      = 250;
+
+# Set the initial state
+x0     = 0;
+
+##############################################################################
+# Generate data
+##############################################################################
+
+(x, y) = helpState.generateData(par, T, x0)
+
+# Plot the measurement
+plt.subplot(3,1,1);
+plt.plot(y,color = '#1B9E77', linewidth=1.5);
+plt.xlabel("time"); plt.ylabel("measurement")
+
+# Plot the states
+plt.subplot(3,1,2);
+plt.plot(x,color = '#D95F02', linewidth=1.5);
+plt.xlabel("time"); plt.ylabel("latent state")
+
+##############################################################################
+# State estimation using the particle filter
+##############################################################################
+
+plt.subplot(3,1,3);
+
+# Using N = 100 particles and plot the estimate of the latent state
+( xhat, ll ) = helpState.pf(y,par,100,T,x0);
+plt.plot(xhat,color = '#7570B3', linewidth=1.5)
+plt.xlabel("time"); plt.ylabel("state estimate")
+
+###################################################################################
+# State estimation using the Kalman filter
+###################################################################################
+
+xhat = helpState.kf(y,par,T,x0,0.01);
+plt.plot(xhat,'k', linewidth=1.5)
+
+##############################################################################
+# End of file
+##############################################################################

python/example2-lgss.py

@@ -0,0 +1,90 @@
+##############################################################################
+#
+# Example of particle Metropolis-Hastings
+# in a linear Gaussian state space model
+#
+# (c) 2015 Johan Dahlin
+# johan.dahlin (at) liu.se
+#
+##############################################################################
+
+# Import some libraries
+import matplotlib.pylab          as plt
+import numpy                     as np
+import stateEstimationHelper     as helpState
+import parameterEstimationHelper as helpParam
+
+# Set the random seed
+np.random.seed( 10 );
+
+##############################################################################
+# Define the model
+##############################################################################
+
+# Here, we use the following model
+#
+# x[tt+1] = phi   * x[tt] + sigmav * v[tt]
+# y[tt]   = x[tt]         + sigmae * e[tt]
+#
+# where v[tt] ~ N(0,1) and e[tt] ~ N(0,1)
+
+# Set the parameters of the model (par[0],par[1],par[2])
+par = np.zeros(3)
+par[0] = 0.75
+par[1] = 1.00
+par[2] = 1.00
+
+# Set the number of time steps to simulate
+T      = 250;
+
+# Set the initial state
+x0     = 0;
+
+##############################################################################
+# Generate data
+##############################################################################
+
+(x, y) = helpState.generateData(par, T, x0)
+
+##############################################################################
+# Parameter estimation using PMH
+##############################################################################
+
+# The inital guess of the parameter
+initPar  = 0.50;
+
+# No. particles in the particle filter ( choose nPart ~ T )
+nPart    = 500;
+
+# The length of the burn-in and the no. iterations of PMH ( nBurnIn < nRuns )
+nBurnIn  = 1000;
+nRuns    = 5000;
+
+# The standard deviation in the random walk proposal
+stepSize = 0.10;
+
+# Run the PMH algorithm
+res = helpParam.pmh(y,initPar,par,nPart,T,x0,helpState.pf,nRuns,stepSize)
+
+
+##############################################################################
+# Plot the results
+##############################################################################
+
+# Plot the parameter posterior estimate
+# Solid black line indicate posterior mean
+plt.subplot(2,1,1);
+plt.hist(res[nBurnIn:nRuns], np.floor(np.sqrt(nRuns-nBurnIn)), normed=1, facecolor='#7570B3')
+plt.xlabel("par[0]"); plt.ylabel("posterior density estimate")
+plt.axvline( np.mean(res[nBurnIn:nRuns]), linewidth=2, color = 'k' )
+
+# Plot the trace of the Markov chain after burn-in
+# Solid black line indicate posterior mean
+plt.subplot(2,1,2);
+plt.plot(np.arange(nBurnIn,nRuns,1),res[nBurnIn:nRuns],color = '#E7298A')
+plt.xlabel("iteration"); plt.ylabel("par[0]")
+plt.axhline( np.mean(res[nBurnIn:nRuns]), linewidth=2, color = 'k' )
+
+##############################################################################
+# End of file
+##############################################################################

python/example3-lgss.py

@@ -0,0 +1,108 @@
+##############################################################################
+#
+# Example of particle Metropolis-Hastings in a stochastic volatility model
+#
+# (c) 2015 Johan Dahlin
+# johan.dahlin (at) liu.se
+#
+##############################################################################
+
+# Import some libraries
+import matplotlib.pylab          as plt
+import numpy                     as np
+import stateEstimationHelper     as helpState
+import parameterEstimationHelper as helpParam
+
+# Set the random seed to replicate results in tutorial
+np.random.seed( 10 );
+
+##############################################################################
+# Define the model
+##############################################################################
+
+# Here, we use the following model
+#
+# x[tt+1] = par[0] * x[tt] + par[1]  * v[tt]
+# y[tt]   = par[2] * exp( xt[tt]/2 ) * e[tt]
+#
+# where v[tt] ~ N(0,1) and e[tt] ~ N(0,1)
+
+# Set the number of time steps
+T      = 500;
+
+##############################################################################
+# Load data
+##############################################################################
+y = np.loadtxt('omxs30data.csv')
+
+##############################################################################
+# Parameter estimation using PMH
+##############################################################################
+
+# The inital guess of the parameter
+initPar  = np.array(( 0.99, 0.13, 0.90));
+
+# No. particles in the particle filter ( choose nPart ~ T )
+nPart    = 500;
+
+# The length of the burn-in and the no. iterations of PMH ( nBurnIn < nRuns )
+nBurnIn  = 2500;
+nRuns    = 7500;
+
+# The standard deviation in the random walk proposal
+stepSize = np.array(( 0.01, 0.05, 0.05));
+
+# Run the PMH algorithm
+res = helpParam.pmh_sv(y,initPar,nPart,T,helpState.pf_sv,nRuns,stepSize)
+
+# State inference
+( xhat, ll ) = helpState.pf_sv(y, np.mean( res[nBurnIn:nRuns,:], axis=0 ) ,nPart,T);
+
+#=============================================================================
+# Plot the results
+#=============================================================================
+
+plt.figure(1);
+
+# Plot the log-returns with volatility estimate
+plt.subplot(4,2,(1,2));
+plt.plot(y,color = '#1B9E77', linewidth=1.5);
+plt.plot(xhat,color = '#D95F02', linewidth=1.5);
+plt.xlabel("time"); plt.ylabel("log-return")
+
+# Plot the parameter posterior estimate
+# Plot the trace of the Markov chain after burn-in
+# Solid black line indicate posterior mean
+plt.subplot(4,2,3);
+plt.hist(res[nBurnIn:nRuns,0], np.floor(np.sqrt(nRuns-nBurnIn)), normed=1, facecolor='#7570B3')
+plt.xlabel("phi"); plt.ylabel("posterior density estimate")
+plt.axvline( np.mean(res[nBurnIn:nRuns,0]), linewidth=1.5, color = 'k' )
+
+plt.subplot(4,2,4);
+plt.plot(np.arange(nBurnIn,nRuns,1),res[nBurnIn:nRuns,0],color='#7570B3')
+plt.xlabel("iteration"); plt.ylabel("trace of phi")
+plt.axhline( np.mean(res[nBurnIn:nRuns,0]), linewidth=1.5, color = 'k' )
+
+plt.subplot(4,2,5);
+plt.hist(res[nBurnIn:nRuns,1], np.floor(np.sqrt(nRuns-nBurnIn)), normed=1, facecolor='#E7298A')
+plt.xlabel("sigmav"); plt.ylabel("posterior density estimate")
+plt.axvline( np.mean(res[nBurnIn:nRuns,1]), linewidth=1.5, color = 'k' )
+
+plt.subplot(4,2,6);
+plt.plot(np.arange(nBurnIn,nRuns,1),res[nBurnIn:nRuns,1],color='#E7298A')
+plt.xlabel("iteration"); plt.ylabel("trace of sigmav")
+plt.axhline( np.mean(res[nBurnIn:nRuns,1]), linewidth=1.5, color = 'k' )
+
+plt.subplot(4,2,7);
+plt.hist(res[nBurnIn:nRuns,2], np.floor(np.sqrt(nRuns-nBurnIn)), normed=1, facecolor='#66A61E')
+plt.xlabel("beta"); plt.ylabel("posterior density estimate")
+plt.axvline( np.mean(res[nBurnIn:nRuns,2]), linewidth=1.5, color = 'k' )
+
+plt.subplot(4,2,8);
+plt.plot(np.arange(nBurnIn,nRuns,1),res[nBurnIn:nRuns,2],color='#66A61E')
+plt.xlabel("iteration"); plt.ylabel("trace of beta")
+plt.axhline( np.mean(res[nBurnIn:nRuns,2]), linewidth=1.5, color = 'k' )
+
+##############################################################################
+# End of file
+##############################################################################

r/example1-lgss.R

@@ -45,7 +45,7 @@ x    = data$x
 y    = data$y
 
 # Export plot to file
-#cairo_pdf("lgss-state.pdf", height = 8, width = 8)
+#cairo_pdf("example1-lgss.pdf", height = 10, width = 8)
 
 # Plot the latent state and observations
 layout(matrix(1:3, 3, 1, byrow = TRUE));

r/example1-lgss.RData.REMOVED.git-id

@@ -0,0 +1 @@
+5f7bb25e089d2065578cf5b15059510548653bc8

r/example2-lgss.R

@@ -69,7 +69,7 @@ res = pmh(data$y,initPar,sigmav,sigmae,nPart,T,x0,nIter,stepSize)
 ##############################################################################
 
 # Export plot to file
-#cairo_pdf("lgss-parameter.pdf", height = 9, width = 8)
+#cairo_pdf("example2-lgss.pdf", height = 10, width = 8)
 
 # Plot the parameter posterior estimate
 # Solid black line indicate posterior mean
@@ -89,6 +89,7 @@ abline(h=mean(res[nBurnIn:nIter]),lwd=1,lty="dotted")
 # Plot the ACF of the Markov chain
 foo=acf( res[nBurnIn:nIter], plot=F)
 plot(foo$lag,foo$acf,col = '#7570B3', type="l",xlab="iteration",ylab="ACF",bty="n",lwd=2)
+polygon(c(foo$lag,rev(foo$lag)),c(foo$acf,rep(0,length(foo$acf))),border=NA,col=rgb(t(col2rgb("#7570B3"))/256,alpha=0.25))
 abline(h=1.96/sqrt(nIter-nBurnIn),lty="dotted")
 abline(h=-1.96/sqrt(nIter-nBurnIn),lty="dotted")
 

r/example2-lgss.RData.REMOVED.git-id

@@ -0,0 +1 @@
+d1941306502b78c61bf5f00ef77ec1cc6403b461