Estimating Parking Occupancy using Smart Meter Transaction Data

Daniel Jordon , Robert C. Hampshire and Tayo Fabusuyi

Abstract.

The excessive search for parking, known as cruising, generates pollution and congestion. Cities are looking for approaches that will reduce the negative impact associated with searching for parking. However, adequately measuring the number of vehicles in search of parking is difficult and requires sensing technologies.

In this paper, we develop an approach that eliminates the need for sensing technology by using parking meter payment transactions to estimate parking occupancy and the number of cars searching for parking. The estimation scheme is based on Particle Markov Chain Monte Carlo. We validate the performance of the Particle Markov Chain Monte Carlo approach using data simulated from a $GI/GI/s$ queue. We show that the approach generates asymptotically unbiased Bayesian estimates of the parking occupancy and underlying model parameters such as arrival rates, average parking time, and the payment compliance rate. Finally, we estimate parking occupancy and cruising using parking meter data from SFpark, a large scale parking experiment and subsequently, compare the Particle Markov Chain Monte Carlo parking occupancy estimates against the ground truth data from the parking sensors. Our approach is easily replicated and scalable given that it only requires using data that cities already possess, namely historical parking payment transactions.

The research was carried out when the authors were all affiliated with the University of Michigan Transportation Research Institute.

1. Introduction

Finding an available parking space is an annoyance of modern life. The search for parking has external costs as well given that excess driving increases traffic congestion and pollutes the environment. In response to these negative externalities, city officials are deploying parking information and pricing systems. Parking occupancy information and market based parking pricing strategies reduce the time to find a parking space and pollution ([fabusuyi2014, pierce2013getting, millard2014curb, chatman2014theory]). These strategies require cities to install both parking occupancy sensors and automated payment stations. However, on-street parking sensors are expensive to install and maintain, and many cities cannot afford to deploy parking sensors. In contrast to parking sensors, the automated payment stations are more affordable and more widely deployed. In this paper, we attempt to answer the following question: is it possible to estimate parking occupancy and the number of drivers searching for parking using only parking meter payment data?

Two challenges stand in the way of answering this question. First, a large fraction of drivers do not pay at all or do not pay enough to cover their parking time. For example, in San Francisco the non-payment rate is roughly 30 percent, driven in part by the perennial problem with the abuse of disabled parking placards for on-street parking [SIRA2014]. Figure 1 plots the probability distribution of the non-payment rates in San Francisco. During metered hours, only 70% of parking time is paid for on the modal block. The second challenge is that drivers almost never pay for the exact amount of parking that they need, as they either pay too much or too little.

Refer to caption — Figure 1. Probability distribution of ratio of paid parking time to occupied time. We observe many blocks experience a low payment compliance rate. Seventy-five percent of the occupied time on the modal block is paid for. Source: San Francisco Metropolitan Transportation Authority

It is these challenges that are estimation approach seeks to address. Specific contributions of our approach include the:

•

The introduction of a modeling framework to infer parking occupancy and the time to find parking given a sequence of parking meter payments using a variant of the Particle Markov Chain Monte Carlo (PMCMC) method. We show that this algorithm generates asymptotically unbiased estimates of parking occupancy of the underlying model parameters, such as arrival rates, average parking time, and non-compliance rates (§4.3).
•

The use of a stochastic queueing model to provide constraints for inferring the unobserved parking occupancy and the model parameters. (§3.1).
•

The validation of the performance of the proposed methods with data from a large scale parking experiment called SFpark (§5.2) and with simulated data from $GI/GI/s$ queue (§5.1).

The sample path construction for the $GI/GI/s$ queue is the foundation of the our queueing inference methods. The sample path formulation leads to a flexible modelling framework for inference. Our proposed method extends to time varying model parameters (arrival rates, service rates, etc). For the ease of exposition, we focus on the case of constant model parameters. Furthermore, our methods extends to the case of random order of service. There exists an analogous sample path construction of a GI/GI/s queue with random order of service (see Sutton and Jordon [sutton2011bayesian]). For the ease of exposition, we focus on the first come first served service discipline.

The balance of the paper is organized as follows: Section 2 provides a review of the existing literature. The details of the application of a stochastic queuing model to parking occupancy is discussed in Section 3. Information on the computational algorithm and how the queuing behavior is inferred from smart parking meter transaction data is provided in Section 4. Numerical experiments are carried out in Section 5 and Section 6, the last section, concludes.

2. Literature Review

Traffic congestion has a huge negative impact not just on local municipalities but also on the global economy. It increases pollution, wastes energy, and costs billions in lost time and productivity. A significant contributor to traffic congestion is vehicles searching for parking and it is estimated that in busy areas, an average of 30 percent of traffic congestion is from drivers cruising for open parking spots [shoup2006cruising]. Also contributing to this problem is the increase in vehicle ownership [Schaller2021], a development which adds to the scarcity of parking spaces, particularly, on-street parking. Monitoring on-street parking is a much more difficult task than monitoring garage parking, which is relatively easy to compute through gate counts of entering and exiting vehicles. And cities often price on-street parking significantly lower than off-street parking, leading drivers rationally to choose to cruise [millard2014curb].

Intelligent transportation systems (ITS) such as smart parking solutions are being used to help relieve such congestion and help drivers navigate this situation by providing real-time and predictive information about current and future traffic and parking availability [fabusuyi2014]. One well-known smart parking solution to the “cruising for parking” problem was SFpark [SIRA2014]. SFpark was a Federal Highway Administration (FHWA) funded parking management project in San Francisco, California from 2009 to 2014. The City of San Francisco installed close to 12,000 magnetometer sensors in 8000 parking spaces, creating a wireless sensor network to collect and distribute information on parking availability [6]. Though conceived as a demonstration project, the problems with replicating the program elsewhere include the enormous cost, which is out of the reach of most cities, the limited life-cycle of the fixed sensors, and the variable pricing mechanism used by SFpark [fabusuyi2018].

Apart from the demonstration program, other approaches of managing on-street parking include sensing using smartphones [nawaz2013] or through sensors attached to vehicles [mathur2010] which requires drivers to opt into the program. Crowdsensing, via on-vehicle sensors, is another opt-in system since it depends entirely on the willingness of drivers to attach the sensors to their vehicles and to provide information for others to use [Liao2016]. Other forms of sensing include a combination of on-vehicle sensors and crowdsourcing [roman2018]. User opt-in is required in this type of project as well.

Detecting parking occupancy is also carried out via vehicle to infrastructure (V2I) communication. An example is a novel on-street parking management system that detects parking occupancy using low-cost Bluetooth beacon transmitters installed in vehicles that are connected to receivers located near on-street parking spots [Chen2019]. Some of these solutions are being advocated by the private sector or by the private sector in partnership with academia. These include the use of radar sensors mounted on streetlights by Siemens called Integrated Smart Parking Solution [Siemens2015]. Here, the radar sensors are mounted on streetlights to scan a bigger area than can be done from a single vehicle. The radar sensors monitor both traffic flows and parking spots, to help provide information to drivers searching for parking

Ford Motor Company and Georgia Institute of Technology (Georgia Tech) have collaborated to build a system that gives information about on-street parking availability using sonars and radars on existing cars [Ford2015]. This technology is available to Ford drivers as a paid service. Other systems use vehicles’ pre-installed parking sensors to better classify on-street areas into either legal or illegal parking spots in order to give more accurate information to drivers searching for parking [coric2013]. Using data from parking meters or kiosks is yet another method of collecting valuable parking information for drivers given that parking meters or kiosks are responsible for 95 percent of on-street paid parking spots [Yang2017]. These transaction data have been used with other data to provide not just the real time information on parking availability but also to forecast short term on-street parking occupancy. Yang et. al. [Yang2019], for example, use a deep neural network -based model that uses parking meter transactions along with traffic and weather conditions to predict parking occupancy. Yang et. al. earlier approach [Yang2017] simulated individual payment and parking behavior using a probabilistic payment model, which when accumulated over all transactions, creates time-varying occupancy estimations with a low error rate.

The use of parking kiosks or meters to estimate parking occupancy is however, complicated by the fact that drivers can neglect to pay for parking, pay more than required for their time spent, or pay less than required. We can use the parking meter payments in a queueing process model of on-street parking to create a queueing inference engine (QIE) [larson1990queue, bertsimas1992deducing, mandelbaum1998estimating]. However, our approach builds on the literature in three significant ways. First, our underlying queueing model has no restriction on the inter-arrival time distribution. Second, our parking meter data makes up a subset of service commencement times, in contrast to the QIE whose observations are the complete set of departure times. Additionally, we only observe a noisy estimate of the service times; namely the payment amounts. Finally, we introduce a completely different framework based on simulating the queuing process and Markov chain Monte Carlo methods.

Our approach is innovative given that the authors know of no existing work that applies a QIE approach for estimating on-street parking spaces using parking meter transaction data. Our approach improves on the San Francisco Municipal Transportation Authority [SIRA2014, demisch2016demand] regression-based model for estimating hourly parking occupancy using payment data. First, our model provides predictions at the time scale of observed payments while SFMTA’s regression-based approach is at the hourly level. Secondly, our framework learns the model parameters on the fly, while previous methods require the user to tune the regression approach for each geographic region. Thus, the ease of replicating our approach is afforded users as they can deploy our method more broadly compared to the regression-based approach.

3. Models

3.1. Parking Model

The unit of analysis is a single on-street parking block with $s$ parking spaces. We model parking occupancy at the block level as a $GI/GI/s$ stochastic queueing model. We define $\{\alpha_{j}\}$ to be a sequence of random independent inter-arrival times of the drivers with a general distribution and $\{\nu_{j}\}$ as a sequence of random independent parking (service) times with a general distribution. These inter-arrival and service times are the random input primitives to the queue. The queue transforms the random input primitives into arrival times, service start times, and departure times. The physics of the queue determine the values of these outputs. The queue transforms the random primitives according to the set of equations (see [kiefer1955theory, krivulin1994recursive, sutton2011bayesian])

	$\displaystyle a_{j}$	$\displaystyle=a_{j-1}+\alpha_{j}$
	$\displaystyle b_{ji}$	$\displaystyle=\max\{d_{j^{\prime}}\|a_{j^{\prime}}<a_{j},p_{j^{\prime}}=i\}$
	$\displaystyle c_{j}$	$\displaystyle=\min_{i\in(0,s]}b_{ji}$
	$\displaystyle p_{j}$	$\displaystyle=\operatorname*{arg\,min}b_{ji}$
	$\displaystyle u_{j}$	$\displaystyle=\max\{a_{j},c_{j}\}$
(1)		$\displaystyle d_{j}$	$\displaystyle=\nu_{j}+u_{j}$

where the outputs are the sequences of arrival times $\{a_{j}\}$ and departure times $\{d_{j}\}$ . Intermediate variables assist the transformation; $u_{j}$ is the service start time of $j$ -th arrival, the $j$ -th arrival parks at space $h_{j}$ , and $c_{j}$ is the first time a parking space becomes available to serve the $j$ -th arrival. This transformation is one-to-one and invertible [sutton2011bayesian]. Additionally, the likelihood of any set of $N$ arrival and departure times is

(2)

p(\mathbf{a},\mathbf{d})=\prod_{j=1}^{N}h_{\alpha}(a_{j}-a_{j-1})\cdot h_{\nu}(d_{j}-u_{j})

where $h_{\alpha}$ is the density of the inter-arrival times with parameters $\theta_{1}$ and $h_{\nu}$ is the density of the service times with parameters $\theta_{2}$ . This product form of the likelihood is due to the independence of the inter-arrival times and the service times.

These equations correspond to drivers parking at available spaces in the order they arrive. We can extend our results to random order of service by slightly modifying the queue input and output equations [sutton2011bayesian]. For ease of presentation, we focus on the first-come-first-served (FCFS) service discipline

Using this transformation, we construct the number of parked cars, the number of cars searching for parking, and the time needed to find an available parking space.

3.2. Payment Observation Model

We now introduce payments into the parking model. We assume that a driver pays immediately after parking. The driver pays at an automated payment station that aggregates payments for all $s$ spaces on the block. After the $k$ -th driver parks, the time remaining on the meter is

m_{k}=\max(m_{k-1}+\beta_{k}-(\tau_{k}-\tau_{k-1}),0)

where $\beta_{k}$ is the amount paid time by the $k$ -th parking driver and $\tau_{k}$ is the time of the $k$ -th payment. The payment amount relates to the true parking time $\nu$ through a general joint distribution. For the sake of simplicity, we assume that the payment amount is the mixture of two distributions. The first component of the mixture is an exponential distribution with a mean equal to the true parking time. The second component of the mixture is a point mass at zero which represents the probability, $p$ , of a driver not paying for parking.

The observations consist of the time of each payment $\tau_{1:K}$ and the time remaining on the meter $m_{1:K}$ . We focus on analyzing a batch of observations $y_{1:K}=(\tau,m)_{1:K}$ .

3.3. State Space

An important feature of the model is that not all parking drivers pay for parking. Therefore, the number of payments do not correspond to the number of arriving vehicles. We now introduce a random variable to track the number of arriving vehicles. Let the random variable $J(k)$ be the number of arriving drivers before the $k$ -th observed payment. Since the $k$ -th payment occurs at time $\tau_{k}$ , the number of arrivals up to the $k$ -th payment satisfies

(3)

\\ J(k)=\operatorname*{arg\,max}_{n}\bigg{\{}\sum_{i=J(k-1)+1}^{n}1_{\{\beta_{i}>0\}}\leq 1\bigg{\}},\quad J(k)=J(k-1)+x_{k}

where $x_{k}$ is a geometric random variable with rate $P(\beta_{k}=0)$

The number of arrivals before the $k$ -th payment, $J(k)$ , evolves according to a Markov chain at payment times. It is a counting process with independent increments. The number of arrivals until $k$ payments has a negative binomial distribution with parameter $1-p$ and $k$ .

The state space is the total number of arrivals and the corresponding collection of arrival times, service start times, and departures times of all customers in the queue at the time of the $k$ -th payment

x_{k}=(j(k),a_{1:j(k)},u_{1:j(k)},d_{1:j(k)}).

The state has two components, the number of arrivals before the $k$ -th payment and the queueing quantities for these arrivals. The dimension of $x_{k}$ grows randomly with $k$ due to the unknown number of arrivals between payments.

4. Inference

Let $\theta$ be the parameters of the inter-arrival time, service time and payment distributions in the parking and payment models. Our goal is to estimate the

(4)

p(\theta,x_{1:K}\,|y_{1:K})=\frac{p_{\theta}(x_{1:K},y_{1:K})\cdot p(\theta)}{p(y_{1:K})}

where $p(\theta)$ is the prior distribution on the model parameters

	$\displaystyle p_{\theta}(x_{1:K},y_{1:K})$	$\displaystyle=p_{\theta}(x_{1:K})\cdot p_{\theta}(y_{1:K}\|x_{1:K})$
		$\displaystyle=\mu_{\theta}{(x_{1})}\prod_{k=2}^{K}f_{\theta}(x_{k}\|x_{k-1})\cdot\prod_{k=1}^{K}g_{\theta}(y_{k}\|x_{k})$

where the process $\{x_{n};n\geq 1\}$ has an initial density $\mu_{\theta}(\cdot)$ and the payments are independent of each other given the state and have a common density $g_{\theta}(y|x)$ where $x$ is the state the system. The transition probability is

f_{\theta}(x_{k}|x_{k-1})=p^{j(k)-j(k-1)-1}(1-p)\cdot\prod_{j=j(k-1)}^{j(k)}h_{\alpha}(a_{j}-a_{j-1})\cdot h_{\nu}(d_{j}-u_{j}).

This joint density $p_{\theta}({x_{1:K},y_{1:K}})$ has the form of a Feynman-Kac model. The rich literature on Feynman-Kac models (see Del Moral [del2000branching]) provides the theoretical foundation for the proposed computational algorithms.

4.1. Particle Filter

The $GI/GI/s$ queue transformation constrains the likely values of the random primitives given the payment observations. First, we consider the case where the parameters of the inter-arrival, service time, and payment distributions are known. Our goal is to compute the posterior distribution of the state given payment observations $p_{\theta}(x_{1:K}|y_{1:K})$ . The traditional Kalman filter [kalman1960new] does not apply to this goal because the dynamics of the state are neither linear nor does it have a Gaussian distribution. The queueing inference engine techniques [larson1990queue, bertsimas1992deducing, mandelbaum1998estimating] do not apply to the goal either due to the general nature of the distribution of the random primitives.

The solution is to simulate sample trajectories from $p_{\theta}(x_{1:K}|y_{1:K})$ using Sequential Monte Carlo (SMC), also known as a particle filter [gordon1993novel, doucet2009tutorial]. The particle filter approximates the posterior distribution over state trajectories using $N$ particles. The particle filter proceeds in three steps. First, for each particle we simulate a one-step state transition. Next, we compute an importance weight for each particle based on the likelihood of observing the current payment given the state of the particle. Finally, we propagate, kill, or replicate each particle in proportion to its importance weight. Particles with larger weights have a higher likelihood of replicating, and those with smaller weights tend to die out. In this way, the surviving particles approximate samples from the posterior distribution $p_{\theta}(x_{1:K}|y_{1:K})$ . For any number of particles, the particle filter produces an unbiased estimate of the posterior distribution [doucet2009tutorial]. It also produces an unbiased estimate of the likelihood of the observations

(5)

\hat{p}(y_{1:K})=\prod_{k=1}^{K}\bigg{(}\frac{1}{N}\sum_{n=1}^{N}w_{k}^{(n)}\bigg{)}

where $w_{k}^{(n)}$ is the importance weight of the $n$ -th particle at step $k$ .

The two ingredients of the particle filter are the state transition probabilities and the observation likelihood model. The particle filter requires only the ability to simulate the one step transitions of the state. The observation likelihood model is more complex to compute than the state transitions.

Data: Payment Observations over time

y_{1:K}

Input:

\theta

, Number of particles

N

, and Effective Sample Size threshold,

N_{ESS}

Sample initial state

\tilde{\mathbf{x}_{0}}\sim p(\mathbf{x}_{0}|\theta)

;

W_{0}^{(n)}\leftarrow\frac{1}{N}

and

L^{0}(\theta)=1

;

for $k=1,\ldots,K$ do

for $n=1,\ldots,N$ do

Sample

\tilde{\mathbf{x}}^{(n)}_{k}\sim f(\mathbf{x}_{k}|\mathbf{x}^{(n)}_{k-1})

using Equation (3.1);

\tilde{\mathbf{x}}^{(n)}_{1:k}\leftarrow(\mathbf{x}^{(n)}_{1:k-1},\tilde{\mathbf{x}}^{(n)}_{k})

;

Sample

\tilde{y}^{(n,h)}_{k}\sim g_{k}(\cdot;\tilde{\mathbf{x}}^{(n)}_{k})

for each

h=1\ldots,H

using the sample path equations;

w_{k}^{(n)}\leftarrow w_{k-1}^{(n)}\cdot g_{\epsilon}(y_{k})

using Equation (6);

L^{k+1}(\theta)\leftarrow L^{k}(\theta)\cdot\frac{1}{N}\sum_{n=1}^{N}w_{k}^{(n)}

;

end for

W_{k}^{(n)}\leftarrow\frac{w_{k}^{(n)}}{\sum_{n=1}^{N}w_{k}^{(n)}}

for each

n=1,\ldots,N

;

if $\sum_{n=1}^{N}\frac{1}{[w_{k}^{(n)}]^{2}}\leq N_{ESS}$ then

Resample particles with replacement according to a multinomial distribution with weights

\{W_{k}\}

. Let

n^{*}

be the resampled index for particle

n

where

P\{n^{*}=l\}=W_{k}^{(n)}

for

l=1,\ldots\,N

. ;

For

n=1,\ldots,N

, set

\mathbf{x}^{(n)}_{1:k}=\tilde{\mathbf{x}}^{n^{*}}_{1:k}

and the weights are uniform

w_{k}^{(n)}=1/N

. ;

else

\mathbf{x}^{(n)}_{1:k}\leftarrow\tilde{\mathbf{x}}^{(n)}_{1:k}

for each

n=1,\ldots,N

;

end if

end for

Output: Distribution of state trajectories over time,

p_{\theta}(\mathbf{x}_{1:k}|y_{1:k})

, and likelihood estimate

\hat{p}_{\theta}(y_{1:K})=\prod_{k=1}^{K}\frac{1}{N}\sum_{n=1}^{N}w_{k}^{(n)}

Algorithm 1 Particle Filter with ABC target.

4.2. Approximate Bayesian Computation (ABC)

In this section, we construct the observation likelihood, $g(y_{k}|\mathbf{x}_{k})$ , of the $k$ -th observation given the state. The observations consist of the payment time and time remaining on the meter just after the $k$ -th payment. The time of the $k$ -th payment must correspond to either an arrival time or a departure time of a previous arrival. However, the structure of the departure times is complex and does not preserve arrival order. Further, the departure times depend on residual service times that are difficult to compute for general service distributions. The amount of the $k$ -th payment depends on the service time of the corresponding driver who need not be the $k$ -th arrival. These factors make the likelihood challenging to compute in closed form, therefore, we simulate an estimate of the likelihood density.

For each particle we simulate $H$ realizations of the $k$ -th payment given the state $x_{k}$ . This simulation is fast due to the sample path construction of the queue and payment models. A particle receives a high weight if the simulated $k$ -th payments are “close” to the observed $k$ -th payment, and receives a low weight otherwise. This is called approximate Bayesian computation (ABC) [marjoram2003markov, beaumont2002approximate]. Jasra et al. [jasra2012filtering] and Calvet et al. [calvet2014accurate] successfully use ABC in the particle filtering context. The full algorithm appears in Algorithm 1. We approximate the likelihood density function by using a kernel-based distance metric in the observation space

(6)

\tilde{g}_{\epsilon}(y_{k})=\frac{1}{{H\cdot\epsilon}}\sum\limits_{h=1}^{H}{K\left({\frac{{{y_{k}}-\tilde{y}^{(h)}_{k}}}{\varepsilon}}\right)}

where $\epsilon$ is the bandwidth of the kernel, a $H$ is the number of simulations.

4.3. Particle Marginal Metropolis-Hastings (PMMH)

We turn to jointly estimating both model parameters. We are interested in the parameters of the inter-arrival time distribution ( $\lambda$ ), the service time distribution ( $\mu$ ) and payment probability ( $p$ ). We implement the Particle Markov Chain Monte Carlo method (PMCMC) [andrieu2010particle]. This combines two powerful methods: Markov Chain Monte Carlo and the particle filter. Specifically, we use the Particle Marginal-Metropolis Hastings (PMMH) method. The outputs of PMMH are an estimate of the data likelihood, $\hat{p}_{\theta}(y_{1:K})$ , and samples from the posterior distribution $p(\theta,x_{1:K}\,|y_{1:K})$ .

The key idea of the approach is that the particle filter generates state proposals within the Metropolis-Hastings algorithm. The particle filter also generates an estimate of the likelihood $p_{\theta}(y_{1:k})$ which is part of the proposal acceptance probability. We accept joint state and parameters proposals that explore parts of the posterior distribution that are more likely. In essence, we are searching over time evolution of the state space and parameters that are the most likely given the observed sequence of payments.

Here is how it works. The input to the method is a batch of payments observations, $y_{1:K}$ . First, we generate proposed values for the parameters using a Gaussian random walk proposal mechanism. Next, we run the particle filter to generate a proposed state trajectory $x_{1:K}$ under the posterior distribution and we generate a likelihood estimate $\hat{p}_{\theta}(y_{1:K})$ under the proposed parameters. We accept the proposed parameters and state trajectory according to an acceptance probability that depends on the likelihood and the parameter proposal distribution. We iterate this process which results in unbiased estimates of joint distribution over the parameters and the state trajectory. The estimates are consistent for any number of particles. The output of PMMH is a distribution of accepted state proposal-parameter pairs. Each accepted proposal of the MCMC step consists of both a state trajectory and corresponding model parameters.

We modify the PMMH algorithm by plugging in the ABC particle filter to generate the state proposals. The full algorithm appears in Algorithm 2. In their discussion of Andrieu et al. [andrieu2010particle], Cornebise and Peters [cornebise2009comments] advocate this approach when the likelihood is intractable. Our main result shows that the ABC-PMMH methods generates unbiased estimates of the the target joint distribution as the number of particles and simulated ABC observations go to infinity.

Theorem 4.1.

The ABC-PMMH converges to the correct posterior distribution when the number of particles grows appropriately with the size of the kernel density estimator, that is

\mathop{\lim}\limits_{N\to\infty}\left|{\log p_{\theta}^{ABC}-\log p_{\theta}}\right|\to 0,

when the number of particles scales with the number of simulated $H$ payments

N=O\left({{{\left({\frac{H}{{\log H}}}\right)}^{-2/5}}}\right)

and the bandwidth of the kernel density estimator decreases with $H$ according to $\varepsilon=c\cdot{\left({\frac{{\log H}}{H}}\right)^{1/5}}$ .

Proof.

The Particle Marginal Metropolis-Hasting acceptance probability depends on the likelihood estimate generated by the particle filter

{p_{\theta}}\left({{y_{1:T}}}\right)=\prod\limits_{t=1}^{T}\bigg{(}\frac{1}{N}\sum\limits_{n=1}^{N}{g_{\theta}\left({y_{t}}|{x_{t}^{n}}\right)\bigg{)}}

where $T$ is the number of steps (payments) and $N$ is the number of particles. This is an unbiased estimator of the likelihood [del2000branching] for any number of particles.

However, the expression for the density function $g$ is intractable. Therefore, we estimate the density $g$ using a Kernel density estimator with bandwidth $\epsilon$ ,

{\hat{g}_{\varepsilon}}\left(y\right)=\frac{1}{{H\cdot\varepsilon}}\sum\limits_{h=1}^{H}{K\left({\frac{{{y_{h}}-y}}{\varepsilon}}\right)}

where $(y_{h})_{h=1}^{H}$ are simulated from true distribution $g$ . The corresponding observation likelihood estimate generated by the particle filter is

\log p_{\theta}^{ABC}=\sum\limits_{t=1}^{T}{\sum\limits_{n=1}^{N}{\frac{1}{{\varepsilon NH}}\sum\limits_{h=1}^{H}{K\left({\tilde{y}_{t}^{n,h};{y_{t}},{\varepsilon^{2}}}\right)}}}.

If $\mathop{\lim}\limits_{N\to\infty}\left|{\log p_{\theta}^{ABC}-\log p_{\theta}}\right|\to 0$ , then the posterior distribution converges to the correct posterior distribution by applying Theorem 4 of Andrieu et al. [andrieu2010particle].

The difference between the true and approximate log likelihood is bounded above

	$\displaystyle\left\|\log p_{\theta}^{ABC}-\log p_{\theta}\right\|$	$\displaystyle\leq\sum_{t=1}^{T}\frac{1}{N}{\sum_{n=1}^{N}{\left\|{\frac{1}{{H\cdot\varepsilon}}\sum_{h=1}^{H}{K\left({\tilde{y}_{k}^{n,h};{y_{k}},{\varepsilon^{2}}}\right)-{g_{\theta}}\left({{y_{k}}\left\|{x_{k}^{n}}\right.}\right)}}\right\|}}$
		$\displaystyle\leq\sum_{t=1}^{T}\frac{1}{N}{\sum_{n=1}^{N}{\mathop{\sup}\limits_{{y_{k}}}\left\|{\frac{1}{{H\cdot\varepsilon}}\sum_{h=1}^{H}{K\left({\tilde{y}_{k}^{n,h};{y_{k}},{\varepsilon^{2}}}\right)-{g_{\theta}}\left({{y_{k}}\left\|{x_{k}^{n}}\right.}\right)}}\right\|}}$
		$\displaystyle\leq T\cdot\delta\cdot{\left({\frac{{\log H}}{H}}\right)^{2/5}}.$

The last inequality follows from the uniform convergence results for kernel density estimators. ∎

Theorem 4.2 (Fan and Yao [fan2003nonlinear] ).

If the density $g$ has a bound second derivative and the bandwidth is chosen so that $\varepsilon=c\cdot{\left({\frac{{\log H}}{H}}\right)^{1/5}}$ , then

\mathop{\sup}\limits_{y}\left|{{{\hat{g}}_{\varepsilon}}\left(y\right)-g\left(y\right)}\right|\leq\delta\cdot{\left({\frac{{\log H}}{H}}\right)^{2/5}}.

We adopt the PMMH method for two reasons. First, the approach allows for flexibility in model specification and only requires the ability to simulate sample paths of the model. Secondly, PMMH handles state spaces with random dimensions. Our state space is of an unknown dimension due to the unknown number of arrivals between payments. In the literature this problem is called trans-dimensional [green2003trans]. One way to perform estimation on trans-dimensional problems is the reversible jump Markov chain Monte Carlo (RJMCMC) method [green1995reversible]. This popular method requires the development of problem specific and often complex proposal moves [karagiannis2013annealed]. There are several papers that use the PMCMC methods in the trans-dimensional problems. For example, PMCMC was used to construct phylogenetic trees [persing2015simulation], and to estimate of volatility in asset pricing models [bauwens2014marginal].

Data: Payment Observations over time

y_{1:K}

Input: Number of particles

N

, Effective Sample Size threshold,

N_{ESS}

, number of Metropolis-Hastings steps

M

for $m=1,\ldots,M$ do

if $m=1\,\,$ then

•

Pick an initial value of parameters: $\mathbf{\theta}(0)$ . ;

•

Run Algorithm 1 targeting $p_{\theta(0)}(\mathbf{x}_{1:K}|y_{1:K})$ and the data likelihood $p_{\theta(0)}(y_{1:K})$ . ;

•

Sample $X_{1:k}(0)\sim\hat{p}_{\theta(0)}(\mathbf{x}_{1:k}|y_{1:K})$ and the likelihood estimate $\hat{p}_{(\theta(0)}(y_{1:K})$ . ;

else

•

Sample from parameter proposals: $\theta^{*}\sim q(\theta|\theta(m-1))$ ;

•

Run Algorithm 1 targeting $p_{\theta^{*}}(\mathbf{x}_{1:j(K}|y_{1:K})$ and the likelihood $p_{\theta^{*}}(y_{1:K})$ . ;

•

Sample $X_{1:k}^{*}\sim\hat{p}_{\theta^{*}}(\cdot|y_{1:K})$ and the likelihood $\hat{p}_{\theta^{*}}(y_{1:K})$ ;

•

Accept proposed parameters and workload trajectories with probability

\min\bigg{(}1,\frac{\hat{p}_{\theta^{*}}(y_{1:K})}{\hat{p}_{\theta(m-1)}(y_{1:K})}\cdot\frac{p(\theta^{*})}{p(\theta(m-1))}\cdot\frac{q(\theta(m-1)|\theta^{*})}{q(\theta^{*}|\theta(m-1))}\bigg{)}.

•

If accepted, set $\theta(m)=\theta^{*},X_{1:k}(m)=X_{1:k}^{*}$ and $\hat{p}_{\theta(m)}(y_{1:K})=\hat{p}_{\theta^{*}}(y_{1:K})$ . If not accepted, set $\theta(m)=\theta(m-1),X_{1:k}(m)=X_{1:k}(m-1)$ and $\hat{p}_{\theta(m)}(y_{1:K})=\hat{p}_{\theta(m-1)}(y_{1:K})$ . ;

end if

end for

Output: Posterior distribution over the parameters and workload trajectories

p(\mathbf{\theta},\mathbf{x}_{1:K}|y_{1:K})

	$\displaystyle\left\|\log p_{\theta}^{ABC}-\log p_{\theta}\right\|$	$\displaystyle\leq\sum_{t=1}^{T}\frac{1}{N}{\sum_{n=1}^{N}{\left\|{\frac{1}{{H\cdot\varepsilon}}\sum_{h=1}^{H}{K\left({\tilde{y}_{k}^{n,h};{y_{k}},{\varepsilon^{2}}}\right)-{g_{\theta}}\left({{y_{k}}\left\|{x_{k}^{n}}\right.}\right)}}\right\|}}$
		$\displaystyle\leq\sum_{t=1}^{T}\frac{1}{N}{\sum_{n=1}^{N}{\mathop{\sup}\limits_{{y_{k}}}\left\|{\frac{1}{{H\cdot\varepsilon}}\sum_{h=1}^{H}{K\left({\tilde{y}_{k}^{n,h};{y_{k}},{\varepsilon^{2}}}\right)-{g_{\theta}}\left({{y_{k}}\left\|{x_{k}^{n}}\right.}\right)}}\right\|}}$
		$\displaystyle\leq T\cdot\delta\cdot{\left({\frac{{\log H}}{H}}\right)^{2/5}}.$

Estimating Parking Occupancy using Smart Meter Transaction Data

Abstract.

1. Introduction

2. Literature Review

3. Models

3.1. Parking Model

3.2. Payment Observation Model

3.3. State Space

4. Inference

4.1. Particle Filter

4.2. Approximate Bayesian Computation (ABC)

4.3. Particle Marginal Metropolis-Hastings (PMMH)

Theorem 4.1.

Proof.

Theorem 4.2 (Fan and Yao [fan2003nonlinear] ).

5. Numerical Experiments

5.1. Simulated Data

5.2. Field Experiment: SFpark

5.3. Sensor Independent Rate Adjustment Approach

5.4. Implications for practice

6. Conclusion

References

Appendix A Joint Distributions