Taking drift-diffusion analysis from the study of turbulent flows to the study of particulate matter smog and air pollutants dynamics
T. Varapongpisan L. Ingsrisawang T. D. Frank


The exact form of the `Ockham factor'
in model selection
    Jonathan Rougier
Jonathan Rougier (corresponding author) is Honorary Professor, School
of Mathematics, University of Bristol (email:
); Carey Priebe is
Professor, Department of Applied Mathematics and Statistics, Johns
Hopkins University (email: ).

       School of Mathematics 

       University of Bristol, UK
       Carey E. Priebe 

       Department of Applied Mathematics and Statistics 

       Johns Hopkins University, USA
     
=====================================================================================================================================================================================================================================================================================================================================================================================================================================================





Drift-diffusion analysis has been introduced in physics 
as a method to study turbulent
flows. 
In the current study, it is proposed to use the method to identify
underlying dynamical models of particulate matter smog, ozone and
nitrogen dioxide concentrations.  
Data from Chiangmai are
considered, which is a major city in the northern part of Thailand that
recently has witnessed a dramatic increase of hospitalization that are assumed
to be related to extreme air pollution levels. 
Three variants of the drift-diffusion analysis method
(kernel-density, binning, linear approximation) are considered. 
It is shown that all three
variants explain the annual pollutant peaks
during the first half of the year by assuming that the parameters of the 
physical-chemical evolution
equations of the pollutants vary periodically throughout the
year. 
Therefore, our analysis provides evidence
that the underlying dynamical models
of the three pollutants being considered are explicitly
time-dependent. 

rift-diffusion analysis, particulate matter, air pollutants 

02.50.Ey, 05.10.Gg, 05.40.-a, 92.60.Sz 



An important task of nonlinear physics and statistics is 
to identify the underlying mechanisms that determine the evolution of systems
on the basis of experimental data.
In this regard, in physics, 
a method has been developed to investigate
turbulent flows <cit.>
that is nowadays frequently called drift-diffusion analysis. The method
was in part motivated by the self-similarity hypothesis of turbulent flows
that in its own merit has been investigated in various systems (see e.g.,
<cit.>).
In recent years,
various studies have 
examined turbulence
using the drift-diffusion analysis 
approach <cit.>,
see also <cit.>.
However, the method turned out to have a broad spectrum of applications (for a
review see <cit.>). 
For example, sport and movement sciences have been taken advantage of
drift-diffusion analysis to identify movement- and posture-related 
dynamical
systems <cit.>. 
Bistable lasers <cit.> and engineering
problems <cit.> have been
examined. 
In what follows, the drift-diffusion analysis approach will be used
to identify underlying laws determining the evolution of air pollutants.
Those laws are assumed to reflect the relevant physical-chemical
evolution equations of the pollutants under consideration as well as
 the impacts of meteorological conditions. 
Monthly extreme value data of air pollutants will be considered
because such extreme air pollutant concentrations 
are assumed to come with serious health 
risks <cit.>
and are likely to increase death  
rates <cit.>. We will analyze data from 
the city of Chiangmai,
Thailand. While Chiangmai is not the largest city of Thailand, it is the
largest city in the northern part of Thailand. Importantly, in recent years,
the number of hospitalizations that are due to high air pollutants 
concentrations is dramatically
increasing in Thailand, in general, and in Chiangmai, in 
particular <cit.>. Therefore, a better understanding of the
dynamics of the monthly extreme scores of air pollutant concentrations would
be beneficial. We will consider the following three air pollutants:
particulate matter that is of 10 micrometers of less (PM_10), ozone
(O_3), and nitrogen dioxide (NO_2). 


Our departure point is a time-discrete sequence of 
observations of pollutant concentrations. This sequence will be 
 referred to as historical trajectory X^h(n) given for the
time points n=1,…,N (with N=60, see below). 
In what follows, n will denote
consecutive months. Our goal
is to derive a stochastic model from the historical trajectory in
analogy to the proposal by Friedrich-Peinke-Renner for historical financial
data and to take seasonal effects into account.
Following the Friedrich-Peinke-Renner method <cit.>, we
consider the increments Y_n(τ) defined by
Y_n(τ) = X(n+τ) - X(n) ⇒ X(n+τ)=Y_n(τ) + X(n)
for τ⩾ 0 and Y_n(0)=0. Parameter τ defines a time scale.
The increments Y_n are assumed to satisfy an evolution equation 
that describes
how
Y_n(τ) evolves
from small scales of a few months (e.g., τ≈ 1,2,3) 
to large scales of a
year (e.g., τ≈ 12,13,14).
In order to determine that evolution equation, 
we consider R increment trajectories of length S with
τ = 0,…,S,  n=1,…,R, and R+S=N. The evolution equation for
Y_n is then obtained using the drift-diffusion 
analysis <cit.>.

Although drift-diffusion 
analysis <cit.>
is as such a non-parametric data analysis method,
it requires to fix a priori the type of the stochastic model
under consideration. In what follows,
we consider a model given in terms of the  
stochastic iterative map 

    Y_n(τ+1)= f(Y_n(τ),m(n,τ)) + g(m(n,τ)) ϵ(τ)  .

In equation (<ref>),
f will be referred to as drift function
(in analogy to the drift function of
a Fokker-Planck equation <cit.>). 
The drift function f is assumed to depend on 
the month m of the year,
where
m  depends on n and τ
like
m(n,τ) = v    v ∈ [1,11]
and 
m(n,τ) = 12    v=0
with v=(n+τ)   12.
In equation (<ref>) 
ϵ(τ) denotes statistically 
independent random variables distributed like a 
normal distribution with mean zero and variance 2. 
The parameter
g⩾ 0 is the noise strength or noise amplitude and, in general, may
depend on the month of the year.
Moreover, g^2 is the noise variance. For the sake of simplicity, it is assumed
that g does not depend on the state Y_n (i.e., an additive noise model
is considered). In early studies by
Friedrich and Peinke <cit.>
and Stanton <cit.>
on the drift-diffusion
analysis, Friedrich, Peinke, and Stanton  
have determined 
representations for drift and diffusion coefficients
of Markov diffusion processes
in terms of conditional averages. In analogy to those
representations, from equation (<ref>) we obtain the
Friedrich-Peinke-Stanton representation of the drift
in terms of the conditional average

    f(z,m) = Y_n(τ+1)|_Y_n(τ)=z .

For the noise variance we obtain

    g^2(m)=1/2[Y_n(τ+1)-f(Y_n(τ),m)]^2  ,

which is not a conditional average because we assume that the noise term 
is
state-independent (i.e., 
additive).
The drift function f can 
approximately be described by means of several methods.
The Friedrich-Peinke binning 
method <cit.>
yields the estimator 

    f(z,m) ≈∑_j=1^K c_j(m) χ_j(z) ,      
    c_j(m) = 
    1/Z_jm∑_τ=0^S [
    ∑_n=1,Y_n(τ) ∈ I_j^R 
    δ_mn Y_n(τ+1)
    ]  ,

where χ_j are indicator functions equal to 1 in appropriately defined
intervals.
We consider Y_n ∈ [A,B) and use K bins of width Δ y
such that y_1=A, y_K+1=B, and y_j=A+(j-1)Δ y. The bin-intervals
are I_j=[y_j,y_j+1). The indicator function is χ_j(z)=1 if 
z ∈ I_j and χ_j(z)=0 otherwise.
In equation (<ref>), δ_mn is the Kronecker function
that equals 1 if the (running) month n corresponds to a particular month m
of the year and zero otherwise. That is, only those
pairs Y_n(τ), Y_n(τ+1) contribute to c_j(m) 
for which the (running) month n is the month of the year m
of interest. Moreover, we have
Z_jm=
∑_τ=0^S [
∑_n=1,Y_n(τ) ∈ I_j^R 
δ_mn].
The kernel density estimation method suggested by 
Stanton <cit.>
yields

    f(z,m) ≈1/Z_m∑_τ=0^S ∑_n=1^R δ_mn Y_n(τ+1)
    exp{ -[z-Y_n(τ)]^2/2h^2}  ,

where the standard deviation h is given by h=s L^-0.2, where
s is the sample standard deviation of all Y_n scores that 
belong to a particular month m of the year 
(i.e., that show up on the sum and for which 
δ_mn=1 holds — 
these are the scores from which the density is estimated) 
and L is the number of such 
scores <cit.>. 
 Moreover, 
Z_m=
∑_τ=0^S ∑_n=1^R δ_mnexp{ -[z-Y_n(τ)]^2/(2h^2) }.
The interpolation modelling method (or regression model method) assumes that
f(z,m) ≈ A_0(m)  + ∑_j=1^p A_j(m) z^j.
For the relative small data sets that will be considered 
below, we will use the model that describes some dependency of f on z 
and features the smallest number of parameters. That is, we 
will consider the order p=1. In this case, equation (<ref>) becomes
the linear regression model 

    f(z,m) ≈ A(m)  + B(m) z  ⇒Y_n(τ+1)= A(m) + B(m) Y_n(τ) + g(m)  ϵ(τ)

with A=A_0 and B=A_1. 
The intercept and slope parameters A(m) and B(m) 
can be estimated 
by fitting the linear regression model equation (<ref>) 
to scatter plots of 
Y_n(τ+1) versus Y_n(τ) given for every month m.
In fact, the 
Y_n(τ+1) versus Y_n(τ) scatter plots are used
to determine f for all three approximations defined by
equations (<ref>), (<ref>) and (<ref>) since 
equations (<ref>), (<ref>) and (<ref>) involve the data pairs
Y_n(τ+1) and Y_n(τ) 
for a fixed month m, that is, all pairs
Y_n(τ+1) and Y_n(τ) for which n corresponds to a particular
month m of the year.
Moreover, from equation (<ref>) it follows that 
g of the linear regression model equation (<ref>) 
can be estimated from the root-mean-squared error RMSE of the
regression model like g(m)=(m)/√(2).






Data were taken from the Pollution Control Department (PCD) of
Thailand <cit.>. Pollutant data 
for 
PM_10, 
O_3, and NO_2 in N=60
months from January 2010 to December 2014 were retrieved for the Provincial
Hall measurement station in Chiangmai. 
Figure <ref> 
shows the pollutant time series. The station measured raw
PM_10 concentrations (in µg/m^3) as averaged values for every day.
From the daily raw data, the PCD determined for each month 
the maximum scores. By contrast,    
O_3 and NO_2, 
raw concentrations (in ppb) were measured by the station every hour. 
From those hourly raw data, maximum scores of the day and 
maximum values for the month were determined. The monthly extreme value
data for 
PM_10, 
O_3, and NO_2 published on the PCD website <cit.> 
were retrieved and analyzed. As mentioned above,
the study of extreme value data is of importance because extreme
pollutant concentrations are related to increased health 
risks <cit.>
and death 
rates <cit.>.
All three pollutants 
PM_10, O_3, and NO_2
showed periodic annual patterns (i.e., seasonal effects), see 
figure <ref>.
PM_10 extreme value concentrations 
peaked in the month of March. Similarly, O_3
extreme value 
concentrations reached maximum values during February, March, and April.
NO_2 extreme value concentrations were the largest in February and March. 
 


For each pollutant trajectory X(n), 
increment trajectories Y_n(τ) were derived 
for reference time points n
in the first three years (i.e.,  n=1,…, R with R=36) such that each
trajectory covered a two years period (i.e., τ=0,…,S with S=23).
From the trajectories Y_n,
scatter plots for each month m showing Y_n(τ+1) versus Y_n(τ) were
obtained. 
From the scatter plots, the drift functions f were determined by means of
 the 3 
different
approximations defined by equations (<ref>), (<ref>) and (<ref>).
Figure <ref> shows the drift functions f(Y,m) thus obtained for
PM_10 for the first four months of the year, January to April. 
The dashed lines
represent diagonals. For January and February, all three approximations of 
f were above the diagonals 
indicating that PM_10 increment concentrations 
increased during those months. That is, if increments were positive in 
January (February),
then they tended to be positive and larger in magnitude in 
February (March). This describes the increase of the PM_10 pollutant
concentration X(n) towards the peaks in March (see figure <ref>A).
By contrast, for March and April, the drift functions
were found to be below the diagonals indicating the PM_10 increment
concentrations
decayed during those months. More precisely, if increments were positive in March (April), 
then they tended to be smaller (closer to zero) or negative in April (May). This corresponds to the decay
of the PM_10 pollutant concentration X(n) from March to May 
[see figure <ref> (A) again]. 





By visual inspection of figure <ref>, 
the kernel density estimation method has the advantage
to account for nonlinear characteristics of f in a smooth fashion. It has
the disadvantage of being described by the whole data sets of 
Y_n(τ+1) and Y_n(τ) pairs 
that contribute to the relevant scatter
plots.  That is, each smooth function f is characterized by a relatively
large set of parameters given in terms of 
Y_n(τ+1) and Y_n(τ) pairs.
The linear approximation has the advantage of being conveniently
characterized by two parameters A and B. It has the disadvantage of not capturing nonlinear effects. The drift function approximation
by means of the binning method can be regarded as a compromise between the two
other approximations. The stair-case like approximate 
drift functions obtained from the
binning method account for nonlinearities. They 
can be described by a bin width
Δ y, the bin centers y_k,c=(y_k+1+y_k)/2, and the function
values f_k. Consequently, to describe f with  K bins, we need
2K+1
parameters. Therefore, the number of parameters to describe f 
with the binning method is larger as
compared to the number of parameters that characterize 
the linear regression model but smaller as compared to the number of 
Y_n(τ+1) and Y_n(τ) parameters 
that are needed to approximate f by 
means of the kernel density estimation method.




The linear regression model approximation, that has the advantage of requiring 
the smallest number of characterizing parameters, was determined in 
detail for all three pollutants. That is, the parameters A, B and g
were estimated as described above. Figure <ref> shows 
the model parameters thus obtained. 
For all three pollutants, the intercept parameter A varied
across the months of the year in a characteristic pattern. 
For PM_10 and O_3 in January and February and 
for NO_2 in January it was positive and assumed the largest positive
values. 
Subsequently, in March and April (PM_10), April and May (O_3 and NO_2),
parameter A was
negative
and assumed the largest-in-the-amount negative values. These patterns, 
as discussed
above in the context of PM_10 
and figure <ref>, describe the mechanism
that leads to the peaks in the original trajectories X(n) around February,
March, and April, see figure <ref>.
For the remaining months from May to December, the
intercept parameters A 
were overall relatively small (i.e., close to zero).
PM_10 and O_3 showed exceptions from this rule in September, where the
parameter values
A
were positive and assumed 20% (PM_10) and 45% (O_3) of their
respective maximal positive A parameters.
For all three pollutants, 
the slope parameter B was found to
be relatively close to unity. 
For PM_10 and NO_2, the noise amplitude g showed  clear 
seasonal peaks around
January, February, and March (PM_10) and March and April (NO_2).






In order to validate the model, we tested the residuals ϵ occurring in 
equation (<ref>).
We determined 
the first
ten lag-s autocorrelation coefficients 
of the residuals  
for each trajectory
Y_n(τ).
The coefficients are shown in 
figure <ref> (panels D, E, F) together with 
single-time-series thresholds <cit.>
(solid lines) for statistically
significant 
and Bonferroni
adjusted <cit.> multiple-tests  thresholds
(dashed lines).
We found that some of the correlation coefficients (in particular, lag-1
correlation coefficients of O_3 and NO_2) 
violated the single-time-series criterion for being not
statistically significant. However, all correlation coefficients were found to
be within the boundaries of the multiple-tests thresholds. Residuals of
trajectories Y_n(τ) were also tested for violation of normality using
the Anderson-Darling normality test. For all PM_10 and
NO_2 trajectories Y_n(τ), the residuals did not violate the normality assumption. For O_3, 
the normality
assumption was violated in 4
out of R=36 trajectories Y_n. 
Overall, the correlation and normality tests 
supported the model 
assumptions.


We showed how to identify stochastic dynamical models for the evolution of air
pollutants on the basis of single, historical trajectories of pollutant
concentrations. To this end, we followed the earlier work on financial data and
considered pollutant increments rather than the raw pollutant data. 
In
addition, three different representation methods of the drift functions of
the dynamical models were used. In doing so, we derived three main results:
First, we found that all three representation methods were consistent with
each other, see figure <ref>. Second,
we were able to show that  experimentally observed 
annual air pollutant peaks were caused by
drift functions of physical-chemical air pollutant systems 
that change qualitatively from the pre-peak
months (e.g., January and February) to the post-peak months (e.g., March and
April), see figure <ref> again. 
These qualitative changes in the drift functions
are assumed to reflect periodic changes in the physical-chemical laws
determining the evolution of the 
PM_10, O_3, and NO_2 pollutant concentrations.
Third, it was found that the linear approximation representation method of 
drift functions (which is the most parsimony method) is sufficient to reproduce
the emergence of the yearly pollutant peaks, see figure <ref>.  

99

friedrich97prl
Friedrich R., Peinke J., Phys. Rev. Lett., 1997, 78, 863, 10.1103/PhysRevLett.78.863.

khomenko10fnl
Khomenko A.V., Lyashenko I.A., Borisyuk V.N., Fluctuation Noise Lett., 
2010, 9, 19,
 10.1142/S0219477510000046.

friedrich97physd
Friedrich R., Peinke J., Physica D, 1997, 102, 147, 10.1016/S0167-2789(96)00235-7.

naert97pre
Naert A., Friedrich R., Peinke J., Phys. Rev. E, 1997, 56, 6719, 10.1103/PhysRevE.56.6719.

tutkun04
Tutkun M., Mydlarski L., New J. Phys., 2004, 6, 49, 10.1088/1367-2630/6/1/049.

kuwahara17
Kuwahara J., Miyata H., Konno H., AIP Conf. Proc., 2017, 1872, 020013, 10.1063/1.4996670.  

friedrich11
Friedrich R., Peinke J., Sahimi M., Tabar M.R.R., Phys. Rep., 2011,
  506, 87, 10.1016/j.physrep.2011.05.003.

davieskurz13
Davies B.L., Kurz M.J., Res. Dev. Disabilities, 2013,
  34, 3648, 10.1016/j.ridd.2013.08.012.

kurz13
Kurz M.J., Arpin D.J., Davies B.L., Harbourne R., Ann. Biomed. Eng., 2013, 41, 1703,
 10.1007/s10439-013-0821-7.

mourik06physletta
Van Mourik A.M., Daffertshofer A., Beek P.J., Phys. Lett. A, 2006,
  351, 13, 10.1016/j.physleta.2005.10.066.

frank06pre
Frank T.D., Friedrich R., Beek P.J., Phys. Rev. E, 2006, 74, 051905, 10.1103/PhysRevE.74.051905.

gottpeinkelippensnagel09pla
Gottschall J., Peinke J., Lippens V., Nagel V., Phys. Lett. A, 2009,
  373, 811, 10.1016/j.physleta.2008.12.026.

frank05jnpcs
Frank T.D., Sondermann M., Ackemann T., Friedrich R., Nonlinear Phenom. Complex
  Syst., 2005, 8, 193.

chiangga10
Chiangga S., Frank T.D., Nonlinear Phenom. Complex Syst., 2010, 13, 32.

gradisek00pre
Gradišek J., Siegert S., Friedrich R., Grabec I., Phys. Rev. E, 2000,
  62, 3146, 10.1103/PhysRevE.62.3146.

lindpeinke14
Lind P.G., Wächter M., Peinke J., J. Phys. Conf. Ser., 2014,
  524, 012179, 10.1088/1742-6596/524/1/012179.

noiray17
Noiray N., J. Eng. Gas Turbines Power, 2016,
  139, 041503, 10.1115/1.4034601.

lid18
Li D., Wang J., Zhang Z., Shen P., Zheng P., Jin M., Lu H., Lin H., Chen K.,
  Environ. Sci. Pollut. Res., 2018, 25, 16135, 10.1007/s11356-018-1759-y.

maji17
Maji K.J., Dikshit A.K., Deshpande A., Environ. Sci. Pollut. Res., 2017, 24, 4709,
 10.1007/s11356-016-8164-1.

pope06
Pope C.A. III, Dockery D.W., J. Air Waste Manage. Assoc., 2006, 56, 709,
 10.1080/10473289.2006.10464485.

giri07
Giri D., Murthy V.K., Adhikary P.R., Khanal S.N., Int. J. Environ. Sci. Technol.,
  2007, 4, 183,
 10.1007/BF03326272.

monterolorenzo11
Montero Lorenzo J.M., Sánchez-Ollero J.L., Fernandez-Aviles G.,
  Int. J. Environ. Res., 2011, 5, 23, 10.22059/ijer.2010.287.

vichit11
Vichit-Vadakan N., Vajanapoom N., Environ. Health Perspect., 2011,
  119, 197, 10.1289/ehp.1103728.

friedrich00aprl
Friedrich R., Peinke J., Renner Ch., Phys. Rev. Lett., 2000, 84, 5224, 10.1103/PhysRevLett.84.5224.

risken89book
Risken H., The Fokker-Planck Equation: Methods of Solution and Applications,
  Springer, Berlin, 1989.

frank05book
Frank T.D., Nonlinear Fokker-Planck Equations: Fundamentals and Applications,
  Springer, Berlin, 2005.

stanton97
Stanton R., J. Finance, 1997, 52, 1973, 10.1111/j.1540-6261.1997.tb02748.x.

silverman86book
Silverman B.W., Density Estimation for Statistics and Data Analysis, Chapman
  and Hall, London, 1986.

frank08aphysa
Frank T.D., Physica A, 2008, 387, 773, 10.1016/j.physa.2007.10.027.


PCD18
PCD (Pollution Control Department, Thailand), Ministry of Natural Resources
  and Environment, Thailand, accessed 2018,
  <http://air4thai.pcd.go.th/webV2/download.php> (in Thai).

diggle90book
Diggle P.J., Time Series: a Biostatistical Introduction, The Clarendon Press,
  Oxford, 1990.

keppel04book
Keppel G., Wickens T.D., Design and Analysis, Pearson Prentice Hall, New
  Jersey, 2004.







	label1
  , ii i, i ,  10900, 
	label2 
	CESPA, ii i,  i, CT 06269, 
	label3 
	ii i,  i, CT 06269, 







-     
	   
	. 
	        
	       ,  
	  .       ׳,       ,      , 
	'     . 
	   - 
	( ,    ).
	,   
	            ,    
	-       .
	,    ,
	   
	  ,   ,      .
	

	
	 ,  ,  
	

