refactorized MATLAB code. Harmonised code and settings between implemenations.

Former-commit-id: 682f903

Author: compops

.gitignore

@@ -1,5 +1,6 @@
 .vscode/
 .R~
+.m~
 
 # Byte-compiled / optimized / DLL files
 __pycache__/

matlab/README.md

@@ -0,0 +1,32 @@
+# MATLAB code for PMH tutorial
+
+This MATLAB code implements the Kalman filter (KF), particle filter (PF) and particle Metropolis-Hastings (PMH) algorithm for two different dynamical models: a linear Gaussian state-space (LGSS) model and a stochastic volatilty (SV) model. Note that the Kalman filter can only be employed for the first of these two models. The details of the code is described in the tutorial paper available at: < http://arxiv.org/pdf/1511.01707 >.
+
+Note that the MATLAB code in this folder covers the basic implementations in the paper. See the R code in r/ for all the implementations and to recreate the results in the tutorial.
+
+Requirements
+--------------
+The code is written and tested for MATLAB 2016b and makes use of the statistics toolbox and the Quandl package. See < https://github.com/quandl/Matlab > for more installation and to download the toolbox. Note that urlread2 is required by the Quandl toolbox and should be installed as detailed in the README file of the Quandl toolbox.
+
+Main script files
+--------------
+These are the main script files that implement the various algorithms discussed in the tutorial.
+
+**example1_lgss.m** State estimation in a LGSS model using the KM and a fully-adapted PF (faPF). The code is discussed in Section 3.1 and the results are presented in Section 3.2 as Figure 4 and Table 1.
+
+**example2_lgss.m** Parameter estimation of one parameter in the LGSS model using PMH with the faPF as the likelihood estimator. The code is discussed in Section 4.1 and the results are presented in Section 4.2 as Figure 5.
+
+**example3_sv.m** Parameter estimation of three parameters in the SV model using PMH with the bootstrap PF as the likelihood estimator. The code is discussed in Section 5.1 and the results are presented in Section 5.2 as Figure 6. The code takes about an hour to run.
+
+Supporting files
+--------------
+**generateData.m** Implements data generation for the LGSS model.
+**kalmanFilter.m** Implements the Kalman filter for the LGSS model.
+**particleFilter.m** Implements the faPF for the LGSS model.
+**particleFilterSVmodel.m** Implements the bPF for the SV model.
+**particleMetropolisHastings.m** Implements the PMH algorithm for the LGSS model.
+**particleMetropolisHastingsSVmodel.m** Implements the PMH algorithm for the SV model.
+
+Adapting the code for another model
+--------------
+See the discussion in *README.MD* in the directory *r/*.

matlab/example1_lgss.m

@@ -1,9 +1,9 @@
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %
-% Example of fully-adapted particle filtering 
-% in a linear Gaussian state space model
+% Example of state estimation in a LGSS model 
+% using particle filters and Kalman filters
 %
-% Copyright (C) 2015 Johan Dahlin < johan.dahlin (at) liu.se >
+% Copyright (C) 2017 Johan Dahlin < liu (at) johandahlin.com.nospam >
 %
 % This program is free software; you can redistribute it and/or modify
 % it under the terms of the GNU General Public License as published by
@@ -21,59 +21,66 @@
 %
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
+% Set random seed
+rng(0)
+
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % 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]
+% x[t + 1] = phi * x[t] + sigmav * v[t]
+% y[t] = x[t] + sigmae * e[t]
 %
-% where v[tt] ~ N(0,1) and e[tt] ~ N(0,1)
+% where v[t] ~ N(0, 1) and e[t] ~ N(0, 1)
 
 % Set the parameters of the model
-phi    = 0.75;
+phi = 0.75;
 sigmav = 1.00;
 sigmae = 0.10;
+theta = [phi sigmav sigmae];
 
 % Set the number of time steps to simulate
-T      = 250;
+T = 250;
 
 % Set the initial state
-x0     = 0;
+initialState = 0;
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % Generate data
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-[x,y] = generateData(phi,sigmav,sigmae,T,x0);
+[x, y] = generateData(theta, T, initialState);
 
 % Plot the measurements and latent states
-subplot(3,1,1); plot(y(2:(T+1)),'LineWidth',1.5,'Color',[27 158 119]/256); 
-xlabel('time'); ylabel('measurement');
+subplot(3,1,1); 
+plot(y(2:(T + 1)), 'LineWidth', 1.5, 'Color', [27 158 119] / 256); 
+xlabel('time'); 
+ylabel('measurement');
 
-subplot(3,1,2); plot(x(2:(T+1)),'LineWidth',1.5,'Color',[217 95 2]/256); 
-xlabel('time'); ylabel('latent state');
+subplot(3,1,2); 
+plot(x(2:(T + 1)), 'LineWidth', 1.5, 'Color', [217 95 2] / 256); 
+xlabel('time'); 
+ylabel('latent state');
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-% State estimation using the particle filter
+% State estimation
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 % Using N = 20 particles and plot the estimate of the latent state
-N      = 20;
-outPF  = sm(y,phi,sigmav,sigmae,N,T,x0);
+outPF  = particleFilter(y, theta, 20, initialState);
 
-subplot(3,1,3); plot(1:(T-1),outPF(2:T),'LineWidth',1.5,'Color',[117 112 179]/256);
-xlabel('time'); ylabel('state estimate');
+% Kalman filter
+outKF  = kalmanFilter(y, theta, initialState, 0.01);
 
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-% State estimation using the Kalman filter
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Plot the difference
+subplot(3,1,3); 
+difference = outPF(2:T) - outKF(2:T);
+plot(1:(T - 1), difference, 'LineWidth', 1.5, 'Color', [117 112 179] / 256);
+xlabel('time'); 
+ylabel('difference in state estimate');
 
-outKF  = kf(y,phi,sigmav,sigmae,T,x0,0.01);
-hold on;
-plot(1:(T-1),outKF(2:T),'k.','LineWidth',1.5);
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % End of file

matlab/example2_lgss.m

@@ -1,10 +1,8 @@
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %
-% Example of PMH in a LGSS model
+% Example of particle Metropolis-Hastings in a LGSS model.
 %
-% Main script
-%
-% Copyright (C) 2015 Johan Dahlin < johan.dahlin (at) liu.se >
+% Copyright (C) 2017 Johan Dahlin < liu (at) johandahlin.com.nospam >
 %
 % This program is free software; you can redistribute it and/or modify
 % it under the terms of the GNU General Public License as published by
@@ -22,86 +20,100 @@
 %
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
+% Set random seed
+rng(0)
+
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % 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]
+% x[t + 1] = phi * x[t] + sigmav * v[t]
+% y[t] = x[t] + sigmae * e[t]
 %
-% where v[tt] ~ N(0,1) and e[tt] ~ N(0,1)
+% where v[t] ~ N(0, 1) and e[t] ~ N(0, 1)
 
 % Set the parameters of the model
-phi    = 0.75;
+phi = 0.75;
 sigmav = 1.00;
 sigmae = 0.10;
+theta = [phi sigmav sigmae];
 
 % Set the number of time steps to simulate
-T      = 500;
+T = 250;
 
 % Set the initial state
-x0     = 0;
+initialState = 0;
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % Generate data
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-[x,y] = generateData(phi,sigmav,sigmae,T,x0);
+[x, y] = generateData(theta, T, initialState);
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % Parameter estimation using PMH
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 % The inital guess of the parameter
-initPar  = 0.50;
+initialPhi = 0.50;
 
-% No. particles in the particle filter ( choose nPart ~ T )
-nPart    = 100;
+% No. particles in the particle filter (choose noParticles ~ T)
+noParticles = 100;
 
-% The length of the burn-in and the no. iterations of PMH ( nBurnIn < nRuns )
-nBurnIn  = 1000;
-nRuns    = 5000;
+% The length of the burn-in and the no. iterations of PMH 
+% (noBurnInIterations < noIterations)
+noBurnInIterations = 1000;
+noIterations = 5000;
 
 % The standard deviation in the random walk proposal
 stepSize = 0.10;
 
 % Run the PMH algorithm
-res = pmh(y,initPar,sigmav,sigmae,nPart,T,x0,nRuns,stepSize);
+res = particleMetropolisHastings(y, initialPhi, [sigmav sigmae],...
+                                 noParticles, initialState,...
+                                 noIterations, stepSize);
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % Plot the results
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
+noBins = floor(sqrt(noIterations - noBurnInIterations));
+grid = noBurnInIterations:noIterations;
+resPhi = res(noBurnInIterations:noIterations);
+
 % Plot the parameter posterior estimate
 % Solid black line indicate posterior mean
-subplot(3,1,1);
-hist(res(nBurnIn:nRuns), floor(sqrt(nRuns-nBurnIn)) );
-xlabel('phi'); ylabel('posterior density estimate');
+subplot(3, 1, 1);
+hist(resPhi, noBins);
+xlabel('phi'); 
+ylabel('posterior density estimate');
 
-h = findobj(gca,'Type','patch');
-set(h,'FaceColor',[117 112 179]/256,'EdgeColor','w');
+h = findobj(gca, 'Type', 'patch');
+set(h, 'FaceColor', [117 112 179] / 256, 'EdgeColor', 'w');
 
 hold on;
-plot([1 1] * mean(res(nBurnIn:nRuns)), [0 500], 'LineWidth', 1);
+plot([1 1] * mean(resPhi), [0 200], 'LineWidth', 3);
 hold off;
 
 % Plot the trace of the Markov chain after burn-in
 % Solid black line indicate posterior mean
-subplot(3,1,2);
-plot(nBurnIn:nRuns,res(nBurnIn:nRuns), 'Color', [117 112 179]/256, 'LineWidth', 1);
-xlabel('iteration'); ylabel('phi');
+subplot(3, 1, 2);
+plot(grid, resPhi, 'Color', [117 112 179] / 256, 'LineWidth', 1);
+xlabel('iteration'); 
+ylabel('phi');
 
 hold on;
-plot([nBurnIn,nRuns], [1 1] * mean(res(nBurnIn:nRuns)), 'k', 'LineWidth', 1);
+plot([grid(1) grid(end)], [1 1] * mean(resPhi), 'k', 'LineWidth', 3);
 hold off;
 
 % Plot ACF of the Markov chain after burn-in
-subplot(3,1,3);
-macf = acf( res(nBurnIn:nRuns), 100, 0 );
-plot(2:101, macf, 'Color', [117 112 179]/256, 'LineWidth', 2);
-xlabel('lag'); ylabel('ACF of phi');
+subplot(3, 1, 3);
+[acf, lags] = xcorr(resPhi - mean(resPhi), 100, 'coeff');
+stem(lags(101:200), acf(101:200), 'Color', [117 112 179] / 256, 'LineWidth', 2);
+xlabel('lag'); 
+ylabel('ACF of phi');
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % End of file