# Bayesian Spatial and Trend Analysis on Ozone Extreme Data in South Korea: 1991–2015

Authors: Cheru Atsmegiorgis Kitabo
Journal: Advances in Meteorology
 (2020)

---

## Abstract

Background. Extreme events like flooding, extreme temperature, and ozone depletion are happening in every corner of the world. Thus, the need to model such rare events having enormous damage has been getting priorities in most countries of the world. Methods. The dataset contains the ozone data from 29 representative air monitoring sites in South Korea collected from 1991 to 2015. Spatial generalized extreme value (GEV) using maximum likelihood estimation (MLE) and two max-stable and Bayesian kriging models are the statistical models used for analysis. Moreover, predictive performances of these statistical models are compared using measures like root-mean-squared error (RMSE), mean absolute error (MAE), relative bias (rBIAS), and relative mean separation (rMSEP) have been utilized. Results. From the time plot of ozone data, extreme ozone concentration is increasing linearly within the specified period. The return level of ozone concentration after 10, 25, 50, and 100 years have been forecasted and showed that there was an increasing trend in ozone extremes. High spatial variability of ozone extreme was observed, and those areas around the territories were having extreme ozone concentration than the centers. Moreover, Bayesian Kriging brought about relatively the minimum RMSE compared to the other models. Conclusion. The extreme ozone concentration has clearly showed a positive trend and spatial variation. Moreover, among the models considered in the paper, the Bayesian Kriging has been chosen as the better model.


---

## Body

## 1. Background

Extreme events like flooding, extreme temperature, and ozone depletion are happening in every corner of the world. Thus, the need to model such rare events is being getting prioritized in most countries of the world. These days, depletion of the stratospheric ozone layer by halogenated chemicals is a global problem. Since the issue was primarily a matter of scientific curiosity during the 1970s and early 1980s, recently, it became an urgent policy question for governments in both developed and developing countries.There are various effects of ozone starting from human health to increasing the chance of existence of natural catastrophe. Some of them are the effects on human and animal health, terrestrial plants, aquatic ecosystems, bio-geochemical cycles, air quality, materials, climate change, and ultraviolet radiation. According to Stedman et al., the effect on human health like eye diseases, skin cancer, and infectious diseases is due to increased penetration of solar UV-B radiation [1]. The changes in plant form and secondary metabolism due to UV-B could have an important implication for plant competitive balance, plant pathogens, and biogeochemical cycles.The climate and environmental applications used multivariate extreme distributions in order to take into account the extremal dependence [2–5]. For the study of extremes in a spatial context, max-stable processes are a relevant framework for modeling and inference. Different parametric models have been later proposed by various authors [6–8]. The extremal t process is a model of that family that extends the Schlather model [9, 10]. It has a supplementary parameter that lets the extremal coefficient function vary in a wider range. Max-stable processes have been used in applications for modeling air pollution using ozone data, rainfall, temperature, snowfall, and snow depth [11–14].Designing the methodology that brings accurate risk measures of these extreme events needed to minimize the casualties on human being and reduce the damage on infrastructure is indispensable. Thus, the importance of modeling such events and getting ready to take precaution measures will have great importance.The main goal of this paper is to model and compare the predictive performance of the extreme values like spatial GEV with MLE estimation, two max-stable spatial models, and the other recent model Bayesian kriging. In most papers, Bayesian kriging has not been considered for comparison. The comparative criterion considered for comparing the predictive performance of Schlather’s, extremalt, and Bayesian kriging is the 5-fold cross validation. The paper considers 29 air monitoring sites found South Korea, and the daily ozone data have been collected from the year 1991 to 2015. 8 hours’ maximum ozone concentration has been selected from the daily maximum ozone concentration. Then, using this daily maximum ozone, the monthly maximum has been selected. From the monthly maximum, annual maximum have been selected. Actually, only ozone season months from May to September have been utilized for extracting data for the model fitting process.


## 2. Methodology

### 2.1. Basics of Extreme Value Theory (EVT)

These days, the statistical modeling of univariate extremes has been well assessed. However, these tools to model spatial extreme data have been active research areas of the time.


#### 2.1.1. Univariate Extreme Value Models

For each station in this study, the univariate extreme value model has been fitted using generalized extreme value (GEV) distribution using Bayesian perspective. The GEV defined as(1)Gz=exp−1+ξz−μσ−1/ξ,on z:1+ξz−μ/σ>0 and with parameters satisfying −∞<μ<∞, σ>0, and −∞<ξ<∞.The maximum value of ozone per siteZt,i∼GEVμt,i,σ,ξ and the parameters of the model are location (μ), scale (σ), and shape (ξ) for each 29 air monitoring sites. Block maxima first approved in work by Fisher and Tippett approach has been utilized [15]. Time series plot of the ozone concentration versus time has been plotted. Trend has been observed. Thus, the location parameter has been expressed as a linear function of time:(2)μt,i=β0,i+β1,it.As the data contain ozone concentration in different sites, spatial dependence is inevitable.


#### 2.1.2. Generalized Extreme Value Distribution

The generalized extreme value (GEV) distribution has a location parameterμ and a scale parameter σ and combines the three types of extreme value distributions by use of a shape parameter ξ. It is usually represented by GEVμ,σ,ξ. The distribution function is as follows: the case ξ<0 corresponds to the Weibull family; for the case ξ=0, the limit of the GEV as ξ⟶0 is used, which leads to the Gumbel family:(3)Gz=exp−exp−z−μσ.Therefore, we now have a single formulation for modeling block maxima data which overcomes the issue of distribution choice from the three types by allowing the data to choose the appropriate tail behavior for the model. The uncertainty for this choice is now incorporated into the model within the uncertainty for the parameter estimate ofξ:(4)Gz=exp−1+ξz−μσ−1/ξ,defined on z:1+ξz−μ/σ>0 and with parameters satisfying −∞<μ<∞, σ>0 and −∞<ξ<∞. The case ξ>0 corresponds to the Fréchet family.The yearly maximum values of the ozone concentration per geographic zonei, Z1,i,…,ZT,i follows generalized extreme value (GEV) distribution of the following form:(5)Zt,i=GEVμt,i,σt,i,ξ,with the location parameter,(6)μt,i=β0,i+β1,it,fort=1,…,25andi=1,2,…,29,where β0,i is the intercept and β1,i is the trend in t.


### 2.2. Return Level Estimation for Block Maxima Approach

Finding the probability of occurrence for large events is the major motivation for modeling extreme data. This demands us to be able to predict the level of the process that we expect to exceed only once in the time period required, which is known as a return level.For yearly maximum data, the value of the return level can be estimated from the GEV distribution function directly. If we letZp be the 1/P year return level of the process, then(7)Gzp=exp−p,where this expression comes from requiring the exceedance rate of zp to average once every 1/P years and considering a constant rate of exceedance. Accordingly, the approximate value of the 1/P year return level for the GEV distribution is given as(8)zp=μ−σξ1−−log1−p−ξ,forξ≠0,μ−σlog−log1−p,forξ=0.


### 2.3. Bayesian Kriging Model Fitting for Annual Maximum Ozone

A spatial model with a simpler version that is considered in this paper utilizes a model with only one latent spatial process and no measurement errors:(9)Level1;Yu=Xβ+σTu,Level2;Tu∼N0,Ry∅,Level3;prβ,σ2,∅.The likelihood function is given by(10)Lβ,σ2,∅|Y∝σ2−n/2Ry∅−1/2exp−12y−Xβ′Ry∅−1y−Xβ.


#### 2.3.1. Posterior for the Mean Parameterβ

The conjugate prior for the mean parameterβ has been considered. Assuming a normal conjugate prior for the mean parameter β, which means(11)β|σ∗2,∅∗∼Nmβ,σ∗2Vβ.The posterior distribution is given by(12)β|Y,σ∗2,∅∗∼NVβ−1+X′σ∗2X−1Vβ−1mβ+X′σ∗2yσ∗2Vβ−1+X′σ∗2X−1∼NβN^;σ∗2VβN^.Thus, normal distribution is a conjugate prior.


#### 2.3.2. Posterior for the Model Parameter forσ2

The posterior distribution forσ2 is obtained by(13)Prσ2|y,β∗,∅∗∝Prσ2Pry|β∗,σ2,∅∗,where the prior distribution is the first term and the second is the likelihood. The posterior distribution for all the three prior distributions is a scaled-inverse-χ2 of the form and is given by(14)σ2|y,β∗,∅∗∼χScI2υ,Q.The choice of prior distribution determines the parametersυ,Q. In a particular case of the inverse-gamma distribution, the scaled-inverse χ2 is the conjugate prior for σ2. This conjugate distribution is specified by two hyperparameters:(15)σ2|β∗,∅∗∼χScI2nσ,Sσ2,which corresponds to nσSσ2/σ2∼χ2nσThe posterior distribution is(16)σ2|y,β∗,∅∗∼χScI2nσ+n,nσSσ2+nσ^2nσ+n,where σ^2=1/ny−Xβ′Ry−1y−Xβ is the maximum likelihood estimator for σ2


#### 2.3.3. Posterior for the Model Parametersβandσ2

The posterior distribution is obtained by(17)Prβ,σ2|y,∅∗∝Prβ,σ2|y,∅∗Pry|β,σ2,∅∗.The posterior distribution is a normal-scaled-inverse-χ2, i.e., the product of normal and scaled-inverse-χ2 densities for both priors considered. It is given as follows:(18)β,σ2|y,ϕ∗∼Nb,VχScI2υ,Q.The prior choice determines the parametersb,v,V,andQ.


### 2.4. Max-Stable Models

A general representation of max-stable processes can be described by two components, a stochastic processXs and a Poisson process Π with the intensity dζ/ζ2 on (0, ∞). Let {Xisi∈N} Xs with EXs=1 and let ζi∈Π be points of the Poisson process. Then, the spectral representation of Schlater’s model is given by(19)Zs=maxi≥1ζiXis,s∈S.This is called Schlater’s spectral representation. A more flexible class of max-stable processes by takingXis to be any stationary Gaussian random field with finite expectation was suggested by Schlather [7]. A stationary max-stable process with unit Fréchet margins can be obtained by(20)Zs=maxiζimax0,Xis,where μ=Emax0,Xis<∞ and ζi denote the points of a Poisson process on (0, ∞) with intensity measure μ−1ζ−1dζ.The correlationρ⋅ and μ−1=2∏ is used if the random process is specified for a stationary isotropic Gaussian random field Xis with unit variance. This process Z(s) is called an extremal Gaussian process, and the bivariate marginal distributions are given by(21)PrZs1≤z1,Zs2≤z2=exp−121z1+1z21+1−2ρh+1z1z2z1+z22,where h is the Euclidean distance between station s1 and s2. For this paper, the power exponential correlation function has been utilized. Extremal dependence was discussed by various scholars in different disciplines [2–4, 16]. Various parametric models have been suggested by different experts [6–8, 17]. The extremal t process is a model of that family that extends the Schlather model [9, 10].


#### 2.4.1. Specification of Parameters of Schlater’s Model

The covariate specification for the location and scale parameters for fitting Schlater’s model is(22)μs=β01+β11latitude+β12longitude,σs=β02+β21latitude2+β22longitude,Shapeξ∼constant.


### 2.5. Model Choice Criteria Using 5-Fold Cross Validation

To compare the precisions of predictions and forecasts obtained from the fitted models, we used some validation criteria [7, 18–20]. We used the root-mean-squared error (RMSE), mean absolute error (MAE), relative BIAS (rBIAS), and relative mean separation (rMSEP). These validation criteria are defined as(23)RMSE=1m∑i=1mz^i−zi2,MAE=1M∑i=1mz^i−zi,rBIAS=1mz¯∑i=1mz^i−zi,rMSEP=∑i=1mz^i−zi2∑i=1mZ¯p−zi2,where m is the total number of observations that we need to validate, zi is the data indexed by i, z^i is the prediction value, and Z¯i and Z¯p are the arithmetic mean of the observations and predictions, respectively.


## 2.1. Basics of Extreme Value Theory (EVT)

These days, the statistical modeling of univariate extremes has been well assessed. However, these tools to model spatial extreme data have been active research areas of the time.


### 2.1.1. Univariate Extreme Value Models

For each station in this study, the univariate extreme value model has been fitted using generalized extreme value (GEV) distribution using Bayesian perspective. The GEV defined as(1)Gz=exp−1+ξz−μσ−1/ξ,on z:1+ξz−μ/σ>0 and with parameters satisfying −∞<μ<∞, σ>0, and −∞<ξ<∞.The maximum value of ozone per siteZt,i∼GEVμt,i,σ,ξ and the parameters of the model are location (μ), scale (σ), and shape (ξ) for each 29 air monitoring sites. Block maxima first approved in work by Fisher and Tippett approach has been utilized [15]. Time series plot of the ozone concentration versus time has been plotted. Trend has been observed. Thus, the location parameter has been expressed as a linear function of time:(2)μt,i=β0,i+β1,it.As the data contain ozone concentration in different sites, spatial dependence is inevitable.


### 2.1.2. Generalized Extreme Value Distribution

The generalized extreme value (GEV) distribution has a location parameterμ and a scale parameter σ and combines the three types of extreme value distributions by use of a shape parameter ξ. It is usually represented by GEVμ,σ,ξ. The distribution function is as follows: the case ξ<0 corresponds to the Weibull family; for the case ξ=0, the limit of the GEV as ξ⟶0 is used, which leads to the Gumbel family:(3)Gz=exp−exp−z−μσ.Therefore, we now have a single formulation for modeling block maxima data which overcomes the issue of distribution choice from the three types by allowing the data to choose the appropriate tail behavior for the model. The uncertainty for this choice is now incorporated into the model within the uncertainty for the parameter estimate ofξ:(4)Gz=exp−1+ξz−μσ−1/ξ,defined on z:1+ξz−μ/σ>0 and with parameters satisfying −∞<μ<∞, σ>0 and −∞<ξ<∞. The case ξ>0 corresponds to the Fréchet family.The yearly maximum values of the ozone concentration per geographic zonei, Z1,i,…,ZT,i follows generalized extreme value (GEV) distribution of the following form:(5)Zt,i=GEVμt,i,σt,i,ξ,with the location parameter,(6)μt,i=β0,i+β1,it,fort=1,…,25andi=1,2,…,29,where β0,i is the intercept and β1,i is the trend in t.


## 2.1.1. Univariate Extreme Value Models

For each station in this study, the univariate extreme value model has been fitted using generalized extreme value (GEV) distribution using Bayesian perspective. The GEV defined as(1)Gz=exp−1+ξz−μσ−1/ξ,on z:1+ξz−μ/σ>0 and with parameters satisfying −∞<μ<∞, σ>0, and −∞<ξ<∞.The maximum value of ozone per siteZt,i∼GEVμt,i,σ,ξ and the parameters of the model are location (μ), scale (σ), and shape (ξ) for each 29 air monitoring sites. Block maxima first approved in work by Fisher and Tippett approach has been utilized [15]. Time series plot of the ozone concentration versus time has been plotted. Trend has been observed. Thus, the location parameter has been expressed as a linear function of time:(2)μt,i=β0,i+β1,it.As the data contain ozone concentration in different sites, spatial dependence is inevitable.


## 2.1.2. Generalized Extreme Value Distribution

The generalized extreme value (GEV) distribution has a location parameterμ and a scale parameter σ and combines the three types of extreme value distributions by use of a shape parameter ξ. It is usually represented by GEVμ,σ,ξ. The distribution function is as follows: the case ξ<0 corresponds to the Weibull family; for the case ξ=0, the limit of the GEV as ξ⟶0 is used, which leads to the Gumbel family:(3)Gz=exp−exp−z−μσ.Therefore, we now have a single formulation for modeling block maxima data which overcomes the issue of distribution choice from the three types by allowing the data to choose the appropriate tail behavior for the model. The uncertainty for this choice is now incorporated into the model within the uncertainty for the parameter estimate ofξ:(4)Gz=exp−1+ξz−μσ−1/ξ,defined on z:1+ξz−μ/σ>0 and with parameters satisfying −∞<μ<∞, σ>0 and −∞<ξ<∞. The case ξ>0 corresponds to the Fréchet family.The yearly maximum values of the ozone concentration per geographic zonei, Z1,i,…,ZT,i follows generalized extreme value (GEV) distribution of the following form:(5)Zt,i=GEVμt,i,σt,i,ξ,with the location parameter,(6)μt,i=β0,i+β1,it,fort=1,…,25andi=1,2,…,29,where β0,i is the intercept and β1,i is the trend in t.


## 2.2. Return Level Estimation for Block Maxima Approach

Finding the probability of occurrence for large events is the major motivation for modeling extreme data. This demands us to be able to predict the level of the process that we expect to exceed only once in the time period required, which is known as a return level.For yearly maximum data, the value of the return level can be estimated from the GEV distribution function directly. If we letZp be the 1/P year return level of the process, then(7)Gzp=exp−p,where this expression comes from requiring the exceedance rate of zp to average once every 1/P years and considering a constant rate of exceedance. Accordingly, the approximate value of the 1/P year return level for the GEV distribution is given as(8)zp=μ−σξ1−−log1−p−ξ,forξ≠0,μ−σlog−log1−p,forξ=0.


## 2.3. Bayesian Kriging Model Fitting for Annual Maximum Ozone

A spatial model with a simpler version that is considered in this paper utilizes a model with only one latent spatial process and no measurement errors:(9)Level1;Yu=Xβ+σTu,Level2;Tu∼N0,Ry∅,Level3;prβ,σ2,∅.The likelihood function is given by(10)Lβ,σ2,∅|Y∝σ2−n/2Ry∅−1/2exp−12y−Xβ′Ry∅−1y−Xβ.


### 2.3.1. Posterior for the Mean Parameterβ

The conjugate prior for the mean parameterβ has been considered. Assuming a normal conjugate prior for the mean parameter β, which means(11)β|σ∗2,∅∗∼Nmβ,σ∗2Vβ.The posterior distribution is given by(12)β|Y,σ∗2,∅∗∼NVβ−1+X′σ∗2X−1Vβ−1mβ+X′σ∗2yσ∗2Vβ−1+X′σ∗2X−1∼NβN^;σ∗2VβN^.Thus, normal distribution is a conjugate prior.


### 2.3.2. Posterior for the Model Parameter forσ2

The posterior distribution forσ2 is obtained by(13)Prσ2|y,β∗,∅∗∝Prσ2Pry|β∗,σ2,∅∗,where the prior distribution is the first term and the second is the likelihood. The posterior distribution for all the three prior distributions is a scaled-inverse-χ2 of the form and is given by(14)σ2|y,β∗,∅∗∼χScI2υ,Q.The choice of prior distribution determines the parametersυ,Q. In a particular case of the inverse-gamma distribution, the scaled-inverse χ2 is the conjugate prior for σ2. This conjugate distribution is specified by two hyperparameters:(15)σ2|β∗,∅∗∼χScI2nσ,Sσ2,which corresponds to nσSσ2/σ2∼χ2nσThe posterior distribution is(16)σ2|y,β∗,∅∗∼χScI2nσ+n,nσSσ2+nσ^2nσ+n,where σ^2=1/ny−Xβ′Ry−1y−Xβ is the maximum likelihood estimator for σ2


### 2.3.3. Posterior for the Model Parametersβandσ2

The posterior distribution is obtained by(17)Prβ,σ2|y,∅∗∝Prβ,σ2|y,∅∗Pry|β,σ2,∅∗.The posterior distribution is a normal-scaled-inverse-χ2, i.e., the product of normal and scaled-inverse-χ2 densities for both priors considered. It is given as follows:(18)β,σ2|y,ϕ∗∼Nb,VχScI2υ,Q.The prior choice determines the parametersb,v,V,andQ.


## 2.3.1. Posterior for the Mean Parameterβ

The conjugate prior for the mean parameterβ has been considered. Assuming a normal conjugate prior for the mean parameter β, which means(11)β|σ∗2,∅∗∼Nmβ,σ∗2Vβ.The posterior distribution is given by(12)β|Y,σ∗2,∅∗∼NVβ−1+X′σ∗2X−1Vβ−1mβ+X′σ∗2yσ∗2Vβ−1+X′σ∗2X−1∼NβN^;σ∗2VβN^.Thus, normal distribution is a conjugate prior.


## 2.3.2. Posterior for the Model Parameter forσ2

The posterior distribution forσ2 is obtained by(13)Prσ2|y,β∗,∅∗∝Prσ2Pry|β∗,σ2,∅∗,where the prior distribution is the first term and the second is the likelihood. The posterior distribution for all the three prior distributions is a scaled-inverse-χ2 of the form and is given by(14)σ2|y,β∗,∅∗∼χScI2υ,Q.The choice of prior distribution determines the parametersυ,Q. In a particular case of the inverse-gamma distribution, the scaled-inverse χ2 is the conjugate prior for σ2. This conjugate distribution is specified by two hyperparameters:(15)σ2|β∗,∅∗∼χScI2nσ,Sσ2,which corresponds to nσSσ2/σ2∼χ2nσThe posterior distribution is(16)σ2|y,β∗,∅∗∼χScI2nσ+n,nσSσ2+nσ^2nσ+n,where σ^2=1/ny−Xβ′Ry−1y−Xβ is the maximum likelihood estimator for σ2


## 2.3.3. Posterior for the Model Parametersβandσ2

The posterior distribution is obtained by(17)Prβ,σ2|y,∅∗∝Prβ,σ2|y,∅∗Pry|β,σ2,∅∗.The posterior distribution is a normal-scaled-inverse-χ2, i.e., the product of normal and scaled-inverse-χ2 densities for both priors considered. It is given as follows:(18)β,σ2|y,ϕ∗∼Nb,VχScI2υ,Q.The prior choice determines the parametersb,v,V,andQ.


## 2.4. Max-Stable Models

A general representation of max-stable processes can be described by two components, a stochastic processXs and a Poisson process Π with the intensity dζ/ζ2 on (0, ∞). Let {Xisi∈N} Xs with EXs=1 and let ζi∈Π be points of the Poisson process. Then, the spectral representation of Schlater’s model is given by(19)Zs=maxi≥1ζiXis,s∈S.This is called Schlater’s spectral representation. A more flexible class of max-stable processes by takingXis to be any stationary Gaussian random field with finite expectation was suggested by Schlather [7]. A stationary max-stable process with unit Fréchet margins can be obtained by(20)Zs=maxiζimax0,Xis,where μ=Emax0,Xis<∞ and ζi denote the points of a Poisson process on (0, ∞) with intensity measure μ−1ζ−1dζ.The correlationρ⋅ and μ−1=2∏ is used if the random process is specified for a stationary isotropic Gaussian random field Xis with unit variance. This process Z(s) is called an extremal Gaussian process, and the bivariate marginal distributions are given by(21)PrZs1≤z1,Zs2≤z2=exp−121z1+1z21+1−2ρh+1z1z2z1+z22,where h is the Euclidean distance between station s1 and s2. For this paper, the power exponential correlation function has been utilized. Extremal dependence was discussed by various scholars in different disciplines [2–4, 16]. Various parametric models have been suggested by different experts [6–8, 17]. The extremal t process is a model of that family that extends the Schlather model [9, 10].


### 2.4.1. Specification of Parameters of Schlater’s Model

The covariate specification for the location and scale parameters for fitting Schlater’s model is(22)μs=β01+β11latitude+β12longitude,σs=β02+β21latitude2+β22longitude,Shapeξ∼constant.


## 2.4.1. Specification of Parameters of Schlater’s Model

The covariate specification for the location and scale parameters for fitting Schlater’s model is(22)μs=β01+β11latitude+β12longitude,σs=β02+β21latitude2+β22longitude,Shapeξ∼constant.


## 2.5. Model Choice Criteria Using 5-Fold Cross Validation

To compare the precisions of predictions and forecasts obtained from the fitted models, we used some validation criteria [7, 18–20]. We used the root-mean-squared error (RMSE), mean absolute error (MAE), relative BIAS (rBIAS), and relative mean separation (rMSEP). These validation criteria are defined as(23)RMSE=1m∑i=1mz^i−zi2,MAE=1M∑i=1mz^i−zi,rBIAS=1mz¯∑i=1mz^i−zi,rMSEP=∑i=1mz^i−zi2∑i=1mZ¯p−zi2,where m is the total number of observations that we need to validate, zi is the data indexed by i, z^i is the prediction value, and Z¯i and Z¯p are the arithmetic mean of the observations and predictions, respectively.


## 3. Results and Discussion

### 3.1. Results

#### 3.1.1. Air Monitoring Sites Selected for the Study

Twenty-five years’ (from 1991 to 2015) extreme ozone data have been collected from the 29 sites. The reason for considering 1991 as the start of data collection was the existence of large number of missing data in almost all 29 sites prior to 1991. The distribution of air monitoring sites in South Korea is displayed in Figure1.Figure 1
The distribution of the selected 29 air monitoring sites in South Korea.


#### 3.1.2. Descriptive Outputs of Maximum Ozone Concentration

The minimum, first quartile, median, third quartile, and maximum are the descriptive measures used to portray the characteristics of extreme ozone concentrations in each 29 sites, and they are depicted in Figure2.Figure 2
Box plot of the distribution of ozone concentration in the 29 air monitoring sites.Relative to the standard ozone concentration in South Korea which is 0.06 ppm, the distribution of maximum ozone concentration in all air monitoring sites was above the standards.Time series plot of the maximum ozone concentration of all stations is displayed in Figure3. As seen in Figure 3, we can clearly observe a positive trend in the ozone concentration with respect to time. This indicates that the location parameter of the generalized extreme value (GEV) distribution can be expressed as a linear combination of time.Figure 3
Time series plots of ozone concentration levels in the 29 air monitoring sites in South Korea.


#### 3.1.3. Station by Station Analysis by GEV

We assumed that the maximum values of O3 per monitoring station, Z1,…,ZT, are independent and follow a GEV distribution of the following form:(24)Zt∼GEVμt,σ,ξ,μt=β0+β1t;t=1,2,…,T=29,where β0 is an intercept parameter while β1 represents a trend in t for the location parameter and t denotes the order in which the measurements were obtained. Moreover, σ and ξ are constants in time and denote the scale and shape parameters of the GEV distribution. The time series plot of all sites in Figure 3 can be taken as a justification to use a linear trend for the location parameter of generalized extreme value distribution of maximum ozone concentration.


#### 3.1.4. Prior Distribution of Parameter of GEV

The prior is defined by the marginal distributionsμ∼N0,102, σ∼logN0,102, and ξ∼N0,102 and the initial value for the prior of trend parameter is set to be 0.728. With the above prior specification, a Markov chain θ0,…,θn with target distribution πθ|x can be generated. Then, chains of length 10,000 have been generated using the initial values θ0=4,1,0.1,0.1 and the proposal standard deviation s=0.02,0.1,0.1,0.1. Convergence of the posteriors of all parameters has been checked using trace plots and Geweke test.As it is depicted in Table1, all the values of β1 for all the stations in the country are positive indicating that there is a positive trend in the maximum ozone concentration with time.Table 1
The slope (β1) and its standard deviation (SD) of the trend line fitted for Bayesian extreme value analysis at 5% level of significance in all 29 air monitoring sites.
Stationsβ1 (SD)Seoul0.014 (0.719)Busan0.011 (0.312)Daegu0.058 (0.787)Incheon0.047 (0.740)Anyang0.017 (0.356)Songnam0.078 (0.810)Uijeongbu0.062 (0.570)Kuangmyeon0.029 (0.450)Ansan0.018 (0.630)Gimcheon0.069 (0.490)Bucheon0.009 (0.060)Chuncehon0.059 (0.740)Kang-neung0.018 (0.390)Wonju0.017 (0.027)Changwon0.038 (0.640)Gwangju0.027 (0.623)Daejon0.083 (0.590)Ulsan0.012 (0.777)Suwon0.010 (0.705)Cheongju0.005 (0.479)Chinju0.027 (0.014)Suncheon0.057 (0.710)Jeonju0.077 (0.640)Kungsan0.007 (0.030)Iksan0.097 (0.840)Yeosu0.128 (0.940)Gwangyang0.076 (0.690)Pohang0.004 (0.150)Gumi0.061 (0.736)


#### 3.1.5. Spatial Generalized Extreme Value Analysis Assuming No Spatial Dependence

In this perspective, the spatial extreme value analysis does not assume spatial dependence among stations in the study areas. Using this specification, the following likelihood is used for this analysis:(25)lZ;ψ=∑i=1n∑j=1klogfGEVZixj,μZj,σZj,ξZj.


#### 3.1.6. Exploratory Plot of Location Parameter versus Geographic Coordinates

We used the exploratory plot in Figure4 to show the relation between the dependent variable ozone extreme value with two covariates latitude and longitude. The plot shows approximately linear relation.Figure 4
The relation between the location parameter of the spatial GEV model and the covariates latitude and longitude.
(a)(b)As seen in the above plot, there seems to be a linear trend between the covariates latitude and longitude and the ozone concentration.


#### 3.1.7. Exploratory Plot of Scale Parameter versus Geographic Coordinates

Similarly, the relationship between the scale parameter of the spatial GEV model and the geographic covariates latitude and longitude has been displayed in Figure5.Figure 5
The relationship between the scale parameter of the spatial GEV model and the covariates latitude and longitude.
(a)(b)


#### 3.1.8. Model Selection Criterion Using Takeuche Information Criteria

Based on the linear trend relation between the ozone concentration with that of latitude and longitude, the spatial generalized extreme value analysis has been fitted. Therefore, a various combinations of covariates have been used to model a spatial GEV model. Using this specification, 9 models have been fitted. Takeuche information criteria (TIC) were used for selection. Based on this criterion, the model having small TIC would be preferable. TIC of the combinations of covariates is displayed in Table2.Table 2
Model selection criteria for spatial GEV model with various combinations of covariates using Takeuche information criteria (TIC).
CombinationsCovariate combinationsTIC values1μ=σ=β01+β11lat+β21lon2738.5852μ=β01+β11lat+β12lon2741.467σ=β02+β21lat2+β22lon3μ=β01+β11lat+β21lon2741.862σ=β02+β1lat4μ=β01+β11lat,σ=β02+β22lon2746.1245μ=β01+β11lat2742.389σ=β02+β21lat+β22lon6μ=β01+β11lon2743.248σ=β02+β21lat+β22lon7μ=β01+β11lat,σ=β02+β12lat2745.9958μ=β01+β11lat,σ=β02+β22lon2745.9959μ=β01+β11lon,σ=β02+β2lon2747.262lon = longitude, lat = latitude, loc = location,μ = location, s = scale, ξ = considered constant.The second combination has been selected as good model for both the location and scale parameters with the consideration of the type of relationships displayed in Figures4 and 5 and actually the smaller TIC value as well. TIC has been used as model selection criteria since Takeuchi [21] proposed Takeuchi information criterion (TIC) as a robust AIC in the spatial GEV model.


#### 3.1.9. Specification of Location and Scale Parameter

The response surface with the covariates latitude and longitude has been fitted as follows:(26)μs=β01+β11lat+β12lon,σs=β02+β21lat2+β22lon,Shapeξ∼constant.Referring the results in Table3, the model coefficient for latitude covariate is negative indicating that there is an inverse relation between ozone concentration and latitude. The negative relation between ozone concentration dispersion and latitude covariate has been observed.Table 3
Spatial GEV response surface coefficients estimation and its standard deviation (SD) using maximum likelihood estimation with covariates latitude and longitude and constant shape parameter.
Location parameterConstant (SD)Latitude (SD)Longitude (SD)−0.051 (0.001)−0.001 (0.00628)0.0013 (0.00185)Scale parameterConstant (SD)Latitude (SD)Shape (SD)0.063 (0.00363)−0.003 (0.00175)−0.4898 (15.6)SD = standard deviation of the parameter.


#### 3.1.10. Interpolation of Return Level for Spatial GEV (No Spatial Dependence)

(1) Interpolation of 5-Year Return Level. For annual maxima data, the value of the return level can be estimated from the GEV distribution function directly.As depicted in Figure6, different colors show the intensity of the extreme ozone concentrations in various air monitoring sites of South Korea. The air monitoring sites found in territories having higher longitude have higher 5-year return level compared to other areas. The 5-year return level is the ozone concentration that is expected to be exceeded once in 5 years.Figure 6
Interpolations of the 5 years’ return level for extreme ozone concentration in 29 air monitoring sites versus longitude (x-axis) and latitude (y-axis) in South Korea.(2) Interpolation of 10 Years’ Return Level. The 10 years’ return levels (Figure 7) have been interpolated for the study areas using the spatial generalized extreme value analysis. A similar pattern has been observed as the 5 years’ return level. Those areas having smaller longitude have a lower return level.Figure 7
Interpolation of the 10 years’ return level for ozone concentration in 29 air monitoring sites versus longitude (x-axis) and latitude (y-axis) in South Korea.The higher return periods have higher return levels. To confirm this relation between return periods and return level, we interpolated the 25, 50, and 100 years’ return levels. These return levels are the magnitude of the extreme ozone concentrations that are expected to be exceeded once in 25, 50, and 100 years, respectively.(3) Interpolation of 25 Years’ Return Level. The areas in the upper left and lower right corner of the 25 years return level plot have higher interpolated return level relative to the other level. Similarly, the areas at the center of the study area South Korea depicted smaller 25 years’ return level. Moreover, areas having longitude between 127.7 and 128.4, irrespective of the latitude, have smaller 25 years’ return level compared to other areas as depicted in Figure 8.Figure 8
Interpolations of the 25 years’ return levels for extreme ozone concentrations in 29 air monitoring sites versus longitude (x-axis) and latitude (y-axis) in South Korea.(4) Interpolation of 50 Years’ Return Level. Similarly, the interpolations of higher return periods like 50 years and 100 years have been analyzed, and accordingly, the results are displayed in Figures 9 and 10, respectively. As seen in the 5, 10, and 25 years’ return levels, the same pattern of interpolated return levels for ozone concentration has been observed for both 50 and 100 years’ return periods.Figure 9
Interpolations of the 50 years’ return levels for extreme ozone concentrations in 29 air monitoring sites versus longitude (x-axis) and latitude (y-axis) in South Korea.Figure 10
Interpolations of the 100 years’ return level for extreme ozone concentrations in 29 air monitoring sites versus longitude (x-axis) and latitude (y-axis) in South Korea.(5) Interpolation of 100 Years’ Return Level


#### 3.1.11. Model Parameter Estimation of Schlater’s Model

The estimated parameter values based on Schlater’s model is displayed in Table4. These values of model parameters are used for interpolation of return levels of extreme ozone in the specified period of time. Moreover, the standard error of each parameter of Schlater’s model depicts the level of dispersion of the possible model parameters in sampling distributions.Table 4
Schlater’s model parameter estimation using the linear correlation between location and scale parameter with that of latitude and longitude (shape parameter considered constant).
Nugget (SD)Range (SD)Smooth (SD)μConstant (SD)μLon (SD)0.003 (0.0012)0.015 (0.0035)61.825 (29.3)−0.049 (0.00053)−0.001 (0.00083)μLat (SD)σCons (SD)σLon (SD)σLat2SDξCons (SD)0.001 (0.0071)−0.064 (0.0061)−0.001 (0.0036)0.001 (0.00083)−0.270 (0.034)Lon = longitude, Lat = latitude, Cons = constant, SD = standard deviation,μ = mean, σ = population standard deviation.


#### 3.1.12. Model Parameter Estimation of Extremalt Model

The estimated model parameter and the corresponding standard error for each parameters value based on the extremalt model specifications have been displayed in Table 5.Table 5
Extremalt max-stable model parameter estimation using the location and scale with geographic coordinates (shape parameter considered constant).
Nugget (SD)Range (SD)Smooth (SD)μConstant (SD)μLon (SD)0.197 (0.00053)0.546 (0.0018)92.897 (18.9)−0.045 (0.00016)−0.001 (0.00073)μLat (SD)σCons (SD)σLon (SD)σLat2SDξCons (SD)0.002 (0.0038)−0.044 (0.013)−0.001 (0.091)0.005 (0.000045)−0.222 (0.026)Lon = longitude, Lat = latitude, Cons = constant, SD = standard deviation,μ= Mean, σ = population standard deviation.


#### 3.1.13. Model Comparison Using 5-Fold Cross Validation

Random classification of these 29 stations in to 5 folds has been done. Four of the folds have 6 air monitoring sites and the last fold has 5 sites. Using 5 years’ return levels as actual ozone concentration and the interpolated ozone concentration as predicted values, the 7 performance measures have been computed for each fold, and we took the average of it (Table6).Table 6
Predictive performance measures of various models using 5 years’ return level as observed maximum ozone value.
Performance measuresSchlather’s modelSpatial GEV (MLE)Bayesian krigingExtremal-t modelMSE0.0060.0140.0040.007RMSE0.0770.1180.0630.083MAE0.0080.0890.0020.009MAPE9.18414.3389.03410.336BIAS−0.008−0.030−0.078−0.009rBIAS−0.092−0.118−0.088−0.104rMSEP1.0041.1821.0031.009We used the root-mean-squared error (RMSE) as a measure of predictive performance. In Table6, using the 5 years’ return level, the RMSE for the Bayesian kriging model brought the minimum RMSE. This indicates that it has good prediction performance relative to the other models. Similar results have been found for the 10, 25, 50, and 100 years’ return levels.


### 3.2. Discussion

Quentin conducted a study with the title “Comparison of spatial extreme value models application to precipitation data” [22]. This paper compares six models from max-stable process and other hierarchical models using the return level. But, the uniqueness of paper mainly lies on the consideration of other means of model comparison called k-fold cross validation in addition to return level measure.When using the 5 years’ return level as an observed value, the RMSE for the Bayesian kriging models brought about the RMSE of 0.063, and it is relatively the minimum compared to 0.077 for Schlater’s model, 0.118 for the spatial GEV with MLE, and 0.083 for the extremalt model. Similarly, using the 25 years’ return level as observed value, the RMSE for Schlater’s, spatial GEV (MLE), Bayesian kriging, and extremal t are 0.083, 0.130, 0.071, and 0.089, respectively. The last return level considered as an observed value was the 50 years’ return level. Using this value, it has brought about an RMSE of 0.089, 0.137, 0.077, and 0.091 for Schlater’s, spatial GEV (MLE), Bayesian kriging, and extremal t models, respectively.Zahidetal conducted a study on spatial prediction and concluded that Bayesian universal kriging fits better than universal kriging [23]. Based on the result from the study by Cui, Bayesian kriging was more precise as compared to those obtained with ordinary kriging [24]. A simulation study highlighted that the methodology of conditional simulations for the extremal t process was effective for a large range of α∈1,6 values which characterizes the spatial asymptotic dependence [25].The finding from the analysis indicated that the Bayesian kriging model brought about the minimum predictive performance measure compared to other extreme value models considered in this research. For researchers interested to extend this research, I would suggest to consider the criterion of amount of time consumed for estimating model parameters and the CPU used for model comparison.


## 3.1. Results

### 3.1.1. Air Monitoring Sites Selected for the Study

Twenty-five years’ (from 1991 to 2015) extreme ozone data have been collected from the 29 sites. The reason for considering 1991 as the start of data collection was the existence of large number of missing data in almost all 29 sites prior to 1991. The distribution of air monitoring sites in South Korea is displayed in Figure1.Figure 1
The distribution of the selected 29 air monitoring sites in South Korea.


### 3.1.2. Descriptive Outputs of Maximum Ozone Concentration

The minimum, first quartile, median, third quartile, and maximum are the descriptive measures used to portray the characteristics of extreme ozone concentrations in each 29 sites, and they are depicted in Figure2.Figure 2
Box plot of the distribution of ozone concentration in the 29 air monitoring sites.Relative to the standard ozone concentration in South Korea which is 0.06 ppm, the distribution of maximum ozone concentration in all air monitoring sites was above the standards.Time series plot of the maximum ozone concentration of all stations is displayed in Figure3. As seen in Figure 3, we can clearly observe a positive trend in the ozone concentration with respect to time. This indicates that the location parameter of the generalized extreme value (GEV) distribution can be expressed as a linear combination of time.Figure 3
Time series plots of ozone concentration levels in the 29 air monitoring sites in South Korea.


### 3.1.3. Station by Station Analysis by GEV

We assumed that the maximum values of O3 per monitoring station, Z1,…,ZT, are independent and follow a GEV distribution of the following form:(24)Zt∼GEVμt,σ,ξ,μt=β0+β1t;t=1,2,…,T=29,where β0 is an intercept parameter while β1 represents a trend in t for the location parameter and t denotes the order in which the measurements were obtained. Moreover, σ and ξ are constants in time and denote the scale and shape parameters of the GEV distribution. The time series plot of all sites in Figure 3 can be taken as a justification to use a linear trend for the location parameter of generalized extreme value distribution of maximum ozone concentration.


### 3.1.4. Prior Distribution of Parameter of GEV

The prior is defined by the marginal distributionsμ∼N0,102, σ∼logN0,102, and ξ∼N0,102 and the initial value for the prior of trend parameter is set to be 0.728. With the above prior specification, a Markov chain θ0,…,θn with target distribution πθ|x can be generated. Then, chains of length 10,000 have been generated using the initial values θ0=4,1,0.1,0.1 and the proposal standard deviation s=0.02,0.1,0.1,0.1. Convergence of the posteriors of all parameters has been checked using trace plots and Geweke test.As it is depicted in Table1, all the values of β1 for all the stations in the country are positive indicating that there is a positive trend in the maximum ozone concentration with time.Table 1
The slope (β1) and its standard deviation (SD) of the trend line fitted for Bayesian extreme value analysis at 5% level of significance in all 29 air monitoring sites.
Stationsβ1 (SD)Seoul0.014 (0.719)Busan0.011 (0.312)Daegu0.058 (0.787)Incheon0.047 (0.740)Anyang0.017 (0.356)Songnam0.078 (0.810)Uijeongbu0.062 (0.570)Kuangmyeon0.029 (0.450)Ansan0.018 (0.630)Gimcheon0.069 (0.490)Bucheon0.009 (0.060)Chuncehon0.059 (0.740)Kang-neung0.018 (0.390)Wonju0.017 (0.027)Changwon0.038 (0.640)Gwangju0.027 (0.623)Daejon0.083 (0.590)Ulsan0.012 (0.777)Suwon0.010 (0.705)Cheongju0.005 (0.479)Chinju0.027 (0.014)Suncheon0.057 (0.710)Jeonju0.077 (0.640)Kungsan0.007 (0.030)Iksan0.097 (0.840)Yeosu0.128 (0.940)Gwangyang0.076 (0.690)Pohang0.004 (0.150)Gumi0.061 (0.736)


### 3.1.5. Spatial Generalized Extreme Value Analysis Assuming No Spatial Dependence

In this perspective, the spatial extreme value analysis does not assume spatial dependence among stations in the study areas. Using this specification, the following likelihood is used for this analysis:(25)lZ;ψ=∑i=1n∑j=1klogfGEVZixj,μZj,σZj,ξZj.


### 3.1.6. Exploratory Plot of Location Parameter versus Geographic Coordinates

We used the exploratory plot in Figure4 to show the relation between the dependent variable ozone extreme value with two covariates latitude and longitude. The plot shows approximately linear relation.Figure 4
The relation between the location parameter of the spatial GEV model and the covariates latitude and longitude.
(a)(b)As seen in the above plot, there seems to be a linear trend between the covariates latitude and longitude and the ozone concentration.


### 3.1.7. Exploratory Plot of Scale Parameter versus Geographic Coordinates

Similarly, the relationship between the scale parameter of the spatial GEV model and the geographic covariates latitude and longitude has been displayed in Figure5.Figure 5
The relationship between the scale parameter of the spatial GEV model and the covariates latitude and longitude.
(a)(b)


### 3.1.8. Model Selection Criterion Using Takeuche Information Criteria

Based on the linear trend relation between the ozone concentration with that of latitude and longitude, the spatial generalized extreme value analysis has been fitted. Therefore, a various combinations of covariates have been used to model a spatial GEV model. Using this specification, 9 models have been fitted. Takeuche information criteria (TIC) were used for selection. Based on this criterion, the model having small TIC would be preferable. TIC of the combinations of covariates is displayed in Table2.Table 2
Model selection criteria for spatial GEV model with various combinations of covariates using Takeuche information criteria (TIC).
CombinationsCovariate combinationsTIC values1μ=σ=β01+β11lat+β21lon2738.5852μ=β01+β11lat+β12lon2741.467σ=β02+β21lat2+β22lon3μ=β01+β11lat+β21lon2741.862σ=β02+β1lat4μ=β01+β11lat,σ=β02+β22lon2746.1245μ=β01+β11lat2742.389σ=β02+β21lat+β22lon6μ=β01+β11lon2743.248σ=β02+β21lat+β22lon7μ=β01+β11lat,σ=β02+β12lat2745.9958μ=β01+β11lat,σ=β02+β22lon2745.9959μ=β01+β11lon,σ=β02+β2lon2747.262lon = longitude, lat = latitude, loc = location,μ = location, s = scale, ξ = considered constant.The second combination has been selected as good model for both the location and scale parameters with the consideration of the type of relationships displayed in Figures4 and 5 and actually the smaller TIC value as well. TIC has been used as model selection criteria since Takeuchi [21] proposed Takeuchi information criterion (TIC) as a robust AIC in the spatial GEV model.


### 3.1.9. Specification of Location and Scale Parameter

The response surface with the covariates latitude and longitude has been fitted as follows:(26)μs=β01+β11lat+β12lon,σs=β02+β21lat2+β22lon,Shapeξ∼constant.Referring the results in Table3, the model coefficient for latitude covariate is negative indicating that there is an inverse relation between ozone concentration and latitude. The negative relation between ozone concentration dispersion and latitude covariate has been observed.Table 3
Spatial GEV response surface coefficients estimation and its standard deviation (SD) using maximum likelihood estimation with covariates latitude and longitude and constant shape parameter.
Location parameterConstant (SD)Latitude (SD)Longitude (SD)−0.051 (0.001)−0.001 (0.00628)0.0013 (0.00185)Scale parameterConstant (SD)Latitude (SD)Shape (SD)0.063 (0.00363)−0.003 (0.00175)−0.4898 (15.6)SD = standard deviation of the parameter.


### 3.1.10. Interpolation of Return Level for Spatial GEV (No Spatial Dependence)

(1) Interpolation of 5-Year Return Level. For annual maxima data, the value of the return level can be estimated from the GEV distribution function directly.As depicted in Figure6, different colors show the intensity of the extreme ozone concentrations in various air monitoring sites of South Korea. The air monitoring sites found in territories having higher longitude have higher 5-year return level compared to other areas. The 5-year return level is the ozone concentration that is expected to be exceeded once in 5 years.Figure 6
Interpolations of the 5 years’ return level for extreme ozone concentration in 29 air monitoring sites versus longitude (x-axis) and latitude (y-axis) in South Korea.(2) Interpolation of 10 Years’ Return Level. The 10 years’ return levels (Figure 7) have been interpolated for the study areas using the spatial generalized extreme value analysis. A similar pattern has been observed as the 5 years’ return level. Those areas having smaller longitude have a lower return level.Figure 7
Interpolation of the 10 years’ return level for ozone concentration in 29 air monitoring sites versus longitude (x-axis) and latitude (y-axis) in South Korea.The higher return periods have higher return levels. To confirm this relation between return periods and return level, we interpolated the 25, 50, and 100 years’ return levels. These return levels are the magnitude of the extreme ozone concentrations that are expected to be exceeded once in 25, 50, and 100 years, respectively.(3) Interpolation of 25 Years’ Return Level. The areas in the upper left and lower right corner of the 25 years return level plot have higher interpolated return level relative to the other level. Similarly, the areas at the center of the study area South Korea depicted smaller 25 years’ return level. Moreover, areas having longitude between 127.7 and 128.4, irrespective of the latitude, have smaller 25 years’ return level compared to other areas as depicted in Figure 8.Figure 8
Interpolations of the 25 years’ return levels for extreme ozone concentrations in 29 air monitoring sites versus longitude (x-axis) and latitude (y-axis) in South Korea.(4) Interpolation of 50 Years’ Return Level. Similarly, the interpolations of higher return periods like 50 years and 100 years have been analyzed, and accordingly, the results are displayed in Figures 9 and 10, respectively. As seen in the 5, 10, and 25 years’ return levels, the same pattern of interpolated return levels for ozone concentration has been observed for both 50 and 100 years’ return periods.Figure 9
Interpolations of the 50 years’ return levels for extreme ozone concentrations in 29 air monitoring sites versus longitude (x-axis) and latitude (y-axis) in South Korea.Figure 10
Interpolations of the 100 years’ return level for extreme ozone concentrations in 29 air monitoring sites versus longitude (x-axis) and latitude (y-axis) in South Korea.(5) Interpolation of 100 Years’ Return Level


### 3.1.11. Model Parameter Estimation of Schlater’s Model

The estimated parameter values based on Schlater’s model is displayed in Table4. These values of model parameters are used for interpolation of return levels of extreme ozone in the specified period of time. Moreover, the standard error of each parameter of Schlater’s model depicts the level of dispersion of the possible model parameters in sampling distributions.Table 4
Schlater’s model parameter estimation using the linear correlation between location and scale parameter with that of latitude and longitude (shape parameter considered constant).
Nugget (SD)Range (SD)Smooth (SD)μConstant (SD)μLon (SD)0.003 (0.0012)0.015 (0.0035)61.825 (29.3)−0.049 (0.00053)−0.001 (0.00083)μLat (SD)σCons (SD)σLon (SD)σLat2SDξCons (SD)0.001 (0.0071)−0.064 (0.0061)−0.001 (0.0036)0.001 (0.00083)−0.270 (0.034)Lon = longitude, Lat = latitude, Cons = constant, SD = standard deviation,μ = mean, σ = population standard deviation.


### 3.1.12. Model Parameter Estimation of Extremalt Model

The estimated model parameter and the corresponding standard error for each parameters value based on the extremalt model specifications have been displayed in Table 5.Table 5
Extremalt max-stable model parameter estimation using the location and scale with geographic coordinates (shape parameter considered constant).
Nugget (SD)Range (SD)Smooth (SD)μConstant (SD)μLon (SD)0.197 (0.00053)0.546 (0.0018)92.897 (18.9)−0.045 (0.00016)−0.001 (0.00073)μLat (SD)σCons (SD)σLon (SD)σLat2SDξCons (SD)0.002 (0.0038)−0.044 (0.013)−0.001 (0.091)0.005 (0.000045)−0.222 (0.026)Lon = longitude, Lat = latitude, Cons = constant, SD = standard deviation,μ= Mean, σ = population standard deviation.


### 3.1.13. Model Comparison Using 5-Fold Cross Validation

Random classification of these 29 stations in to 5 folds has been done. Four of the folds have 6 air monitoring sites and the last fold has 5 sites. Using 5 years’ return levels as actual ozone concentration and the interpolated ozone concentration as predicted values, the 7 performance measures have been computed for each fold, and we took the average of it (Table6).Table 6
Predictive performance measures of various models using 5 years’ return level as observed maximum ozone value.
Performance measuresSchlather’s modelSpatial GEV (MLE)Bayesian krigingExtremal-t modelMSE0.0060.0140.0040.007RMSE0.0770.1180.0630.083MAE0.0080.0890.0020.009MAPE9.18414.3389.03410.336BIAS−0.008−0.030−0.078−0.009rBIAS−0.092−0.118−0.088−0.104rMSEP1.0041.1821.0031.009We used the root-mean-squared error (RMSE) as a measure of predictive performance. In Table6, using the 5 years’ return level, the RMSE for the Bayesian kriging model brought the minimum RMSE. This indicates that it has good prediction performance relative to the other models. Similar results have been found for the 10, 25, 50, and 100 years’ return levels.


## 3.1.1. Air Monitoring Sites Selected for the Study

Twenty-five years’ (from 1991 to 2015) extreme ozone data have been collected from the 29 sites. The reason for considering 1991 as the start of data collection was the existence of large number of missing data in almost all 29 sites prior to 1991. The distribution of air monitoring sites in South Korea is displayed in Figure1.Figure 1
The distribution of the selected 29 air monitoring sites in South Korea.


## 3.1.2. Descriptive Outputs of Maximum Ozone Concentration

The minimum, first quartile, median, third quartile, and maximum are the descriptive measures used to portray the characteristics of extreme ozone concentrations in each 29 sites, and they are depicted in Figure2.Figure 2
Box plot of the distribution of ozone concentration in the 29 air monitoring sites.Relative to the standard ozone concentration in South Korea which is 0.06 ppm, the distribution of maximum ozone concentration in all air monitoring sites was above the standards.Time series plot of the maximum ozone concentration of all stations is displayed in Figure3. As seen in Figure 3, we can clearly observe a positive trend in the ozone concentration with respect to time. This indicates that the location parameter of the generalized extreme value (GEV) distribution can be expressed as a linear combination of time.Figure 3
Time series plots of ozone concentration levels in the 29 air monitoring sites in South Korea.


## 3.1.3. Station by Station Analysis by GEV

We assumed that the maximum values of O3 per monitoring station, Z1,…,ZT, are independent and follow a GEV distribution of the following form:(24)Zt∼GEVμt,σ,ξ,μt=β0+β1t;t=1,2,…,T=29,where β0 is an intercept parameter while β1 represents a trend in t for the location parameter and t denotes the order in which the measurements were obtained. Moreover, σ and ξ are constants in time and denote the scale and shape parameters of the GEV distribution. The time series plot of all sites in Figure 3 can be taken as a justification to use a linear trend for the location parameter of generalized extreme value distribution of maximum ozone concentration.


## 3.1.4. Prior Distribution of Parameter of GEV

The prior is defined by the marginal distributionsμ∼N0,102, σ∼logN0,102, and ξ∼N0,102 and the initial value for the prior of trend parameter is set to be 0.728. With the above prior specification, a Markov chain θ0,…,θn with target distribution πθ|x can be generated. Then, chains of length 10,000 have been generated using the initial values θ0=4,1,0.1,0.1 and the proposal standard deviation s=0.02,0.1,0.1,0.1. Convergence of the posteriors of all parameters has been checked using trace plots and Geweke test.As it is depicted in Table1, all the values of β1 for all the stations in the country are positive indicating that there is a positive trend in the maximum ozone concentration with time.Table 1
The slope (β1) and its standard deviation (SD) of the trend line fitted for Bayesian extreme value analysis at 5% level of significance in all 29 air monitoring sites.
Stationsβ1 (SD)Seoul0.014 (0.719)Busan0.011 (0.312)Daegu0.058 (0.787)Incheon0.047 (0.740)Anyang0.017 (0.356)Songnam0.078 (0.810)Uijeongbu0.062 (0.570)Kuangmyeon0.029 (0.450)Ansan0.018 (0.630)Gimcheon0.069 (0.490)Bucheon0.009 (0.060)Chuncehon0.059 (0.740)Kang-neung0.018 (0.390)Wonju0.017 (0.027)Changwon0.038 (0.640)Gwangju0.027 (0.623)Daejon0.083 (0.590)Ulsan0.012 (0.777)Suwon0.010 (0.705)Cheongju0.005 (0.479)Chinju0.027 (0.014)Suncheon0.057 (0.710)Jeonju0.077 (0.640)Kungsan0.007 (0.030)Iksan0.097 (0.840)Yeosu0.128 (0.940)Gwangyang0.076 (0.690)Pohang0.004 (0.150)Gumi0.061 (0.736)


## 3.1.5. Spatial Generalized Extreme Value Analysis Assuming No Spatial Dependence

In this perspective, the spatial extreme value analysis does not assume spatial dependence among stations in the study areas. Using this specification, the following likelihood is used for this analysis:(25)lZ;ψ=∑i=1n∑j=1klogfGEVZixj,μZj,σZj,ξZj.


## 3.1.6. Exploratory Plot of Location Parameter versus Geographic Coordinates

We used the exploratory plot in Figure4 to show the relation between the dependent variable ozone extreme value with two covariates latitude and longitude. The plot shows approximately linear relation.Figure 4
The relation between the location parameter of the spatial GEV model and the covariates latitude and longitude.
(a)(b)As seen in the above plot, there seems to be a linear trend between the covariates latitude and longitude and the ozone concentration.


## 3.1.7. Exploratory Plot of Scale Parameter versus Geographic Coordinates

Similarly, the relationship between the scale parameter of the spatial GEV model and the geographic covariates latitude and longitude has been displayed in Figure5.Figure 5
The relationship between the scale parameter of the spatial GEV model and the covariates latitude and longitude.
(a)(b)


## 3.1.8. Model Selection Criterion Using Takeuche Information Criteria

Based on the linear trend relation between the ozone concentration with that of latitude and longitude, the spatial generalized extreme value analysis has been fitted. Therefore, a various combinations of covariates have been used to model a spatial GEV model. Using this specification, 9 models have been fitted. Takeuche information criteria (TIC) were used for selection. Based on this criterion, the model having small TIC would be preferable. TIC of the combinations of covariates is displayed in Table2.Table 2
Model selection criteria for spatial GEV model with various combinations of covariates using Takeuche information criteria (TIC).
CombinationsCovariate combinationsTIC values1μ=σ=β01+β11lat+β21lon2738.5852μ=β01+β11lat+β12lon2741.467σ=β02+β21lat2+β22lon3μ=β01+β11lat+β21lon2741.862σ=β02+β1lat4μ=β01+β11lat,σ=β02+β22lon2746.1245μ=β01+β11lat2742.389σ=β02+β21lat+β22lon6μ=β01+β11lon2743.248σ=β02+β21lat+β22lon7μ=β01+β11lat,σ=β02+β12lat2745.9958μ=β01+β11lat,σ=β02+β22lon2745.9959μ=β01+β11lon,σ=β02+β2lon2747.262lon = longitude, lat = latitude, loc = location,μ = location, s = scale, ξ = considered constant.The second combination has been selected as good model for both the location and scale parameters with the consideration of the type of relationships displayed in Figures4 and 5 and actually the smaller TIC value as well. TIC has been used as model selection criteria since Takeuchi [21] proposed Takeuchi information criterion (TIC) as a robust AIC in the spatial GEV model.


## 3.1.9. Specification of Location and Scale Parameter

The response surface with the covariates latitude and longitude has been fitted as follows:(26)μs=β01+β11lat+β12lon,σs=β02+β21lat2+β22lon,Shapeξ∼constant.Referring the results in Table3, the model coefficient for latitude covariate is negative indicating that there is an inverse relation between ozone concentration and latitude. The negative relation between ozone concentration dispersion and latitude covariate has been observed.Table 3
Spatial GEV response surface coefficients estimation and its standard deviation (SD) using maximum likelihood estimation with covariates latitude and longitude and constant shape parameter.
Location parameterConstant (SD)Latitude (SD)Longitude (SD)−0.051 (0.001)−0.001 (0.00628)0.0013 (0.00185)Scale parameterConstant (SD)Latitude (SD)Shape (SD)0.063 (0.00363)−0.003 (0.00175)−0.4898 (15.6)SD = standard deviation of the parameter.


## 3.1.10. Interpolation of Return Level for Spatial GEV (No Spatial Dependence)

(1) Interpolation of 5-Year Return Level. For annual maxima data, the value of the return level can be estimated from the GEV distribution function directly.As depicted in Figure6, different colors show the intensity of the extreme ozone concentrations in various air monitoring sites of South Korea. The air monitoring sites found in territories having higher longitude have higher 5-year return level compared to other areas. The 5-year return level is the ozone concentration that is expected to be exceeded once in 5 years.Figure 6
Interpolations of the 5 years’ return level for extreme ozone concentration in 29 air monitoring sites versus longitude (x-axis) and latitude (y-axis) in South Korea.(2) Interpolation of 10 Years’ Return Level. The 10 years’ return levels (Figure 7) have been interpolated for the study areas using the spatial generalized extreme value analysis. A similar pattern has been observed as the 5 years’ return level. Those areas having smaller longitude have a lower return level.Figure 7
Interpolation of the 10 years’ return level for ozone concentration in 29 air monitoring sites versus longitude (x-axis) and latitude (y-axis) in South Korea.The higher return periods have higher return levels. To confirm this relation between return periods and return level, we interpolated the 25, 50, and 100 years’ return levels. These return levels are the magnitude of the extreme ozone concentrations that are expected to be exceeded once in 25, 50, and 100 years, respectively.(3) Interpolation of 25 Years’ Return Level. The areas in the upper left and lower right corner of the 25 years return level plot have higher interpolated return level relative to the other level. Similarly, the areas at the center of the study area South Korea depicted smaller 25 years’ return level. Moreover, areas having longitude between 127.7 and 128.4, irrespective of the latitude, have smaller 25 years’ return level compared to other areas as depicted in Figure 8.Figure 8
Interpolations of the 25 years’ return levels for extreme ozone concentrations in 29 air monitoring sites versus longitude (x-axis) and latitude (y-axis) in South Korea.(4) Interpolation of 50 Years’ Return Level. Similarly, the interpolations of higher return periods like 50 years and 100 years have been analyzed, and accordingly, the results are displayed in Figures 9 and 10, respectively. As seen in the 5, 10, and 25 years’ return levels, the same pattern of interpolated return levels for ozone concentration has been observed for both 50 and 100 years’ return periods.Figure 9
Interpolations of the 50 years’ return levels for extreme ozone concentrations in 29 air monitoring sites versus longitude (x-axis) and latitude (y-axis) in South Korea.Figure 10
Interpolations of the 100 years’ return level for extreme ozone concentrations in 29 air monitoring sites versus longitude (x-axis) and latitude (y-axis) in South Korea.(5) Interpolation of 100 Years’ Return Level


## 3.1.11. Model Parameter Estimation of Schlater’s Model

The estimated parameter values based on Schlater’s model is displayed in Table4. These values of model parameters are used for interpolation of return levels of extreme ozone in the specified period of time. Moreover, the standard error of each parameter of Schlater’s model depicts the level of dispersion of the possible model parameters in sampling distributions.Table 4
Schlater’s model parameter estimation using the linear correlation between location and scale parameter with that of latitude and longitude (shape parameter considered constant).
Nugget (SD)Range (SD)Smooth (SD)μConstant (SD)μLon (SD)0.003 (0.0012)0.015 (0.0035)61.825 (29.3)−0.049 (0.00053)−0.001 (0.00083)μLat (SD)σCons (SD)σLon (SD)σLat2SDξCons (SD)0.001 (0.0071)−0.064 (0.0061)−0.001 (0.0036)0.001 (0.00083)−0.270 (0.034)Lon = longitude, Lat = latitude, Cons = constant, SD = standard deviation,μ = mean, σ = population standard deviation.


## 3.1.12. Model Parameter Estimation of Extremalt Model

The estimated model parameter and the corresponding standard error for each parameters value based on the extremalt model specifications have been displayed in Table 5.Table 5
Extremalt max-stable model parameter estimation using the location and scale with geographic coordinates (shape parameter considered constant).
Nugget (SD)Range (SD)Smooth (SD)μConstant (SD)μLon (SD)0.197 (0.00053)0.546 (0.0018)92.897 (18.9)−0.045 (0.00016)−0.001 (0.00073)μLat (SD)σCons (SD)σLon (SD)σLat2SDξCons (SD)0.002 (0.0038)−0.044 (0.013)−0.001 (0.091)0.005 (0.000045)−0.222 (0.026)Lon = longitude, Lat = latitude, Cons = constant, SD = standard deviation,μ= Mean, σ = population standard deviation.


## 3.1.13. Model Comparison Using 5-Fold Cross Validation

Random classification of these 29 stations in to 5 folds has been done. Four of the folds have 6 air monitoring sites and the last fold has 5 sites. Using 5 years’ return levels as actual ozone concentration and the interpolated ozone concentration as predicted values, the 7 performance measures have been computed for each fold, and we took the average of it (Table6).Table 6
Predictive performance measures of various models using 5 years’ return level as observed maximum ozone value.
Performance measuresSchlather’s modelSpatial GEV (MLE)Bayesian krigingExtremal-t modelMSE0.0060.0140.0040.007RMSE0.0770.1180.0630.083MAE0.0080.0890.0020.009MAPE9.18414.3389.03410.336BIAS−0.008−0.030−0.078−0.009rBIAS−0.092−0.118−0.088−0.104rMSEP1.0041.1821.0031.009We used the root-mean-squared error (RMSE) as a measure of predictive performance. In Table6, using the 5 years’ return level, the RMSE for the Bayesian kriging model brought the minimum RMSE. This indicates that it has good prediction performance relative to the other models. Similar results have been found for the 10, 25, 50, and 100 years’ return levels.


## 3.2. Discussion

Quentin conducted a study with the title “Comparison of spatial extreme value models application to precipitation data” [22]. This paper compares six models from max-stable process and other hierarchical models using the return level. But, the uniqueness of paper mainly lies on the consideration of other means of model comparison called k-fold cross validation in addition to return level measure.When using the 5 years’ return level as an observed value, the RMSE for the Bayesian kriging models brought about the RMSE of 0.063, and it is relatively the minimum compared to 0.077 for Schlater’s model, 0.118 for the spatial GEV with MLE, and 0.083 for the extremalt model. Similarly, using the 25 years’ return level as observed value, the RMSE for Schlater’s, spatial GEV (MLE), Bayesian kriging, and extremal t are 0.083, 0.130, 0.071, and 0.089, respectively. The last return level considered as an observed value was the 50 years’ return level. Using this value, it has brought about an RMSE of 0.089, 0.137, 0.077, and 0.091 for Schlater’s, spatial GEV (MLE), Bayesian kriging, and extremal t models, respectively.Zahidetal conducted a study on spatial prediction and concluded that Bayesian universal kriging fits better than universal kriging [23]. Based on the result from the study by Cui, Bayesian kriging was more precise as compared to those obtained with ordinary kriging [24]. A simulation study highlighted that the methodology of conditional simulations for the extremal t process was effective for a large range of α∈1,6 values which characterizes the spatial asymptotic dependence [25].The finding from the analysis indicated that the Bayesian kriging model brought about the minimum predictive performance measure compared to other extreme value models considered in this research. For researchers interested to extend this research, I would suggest to consider the criterion of amount of time consumed for estimating model parameters and the CPU used for model comparison.


## 4. Conclusion

In this paper, spatial GEV using MLE estimation, two models from max-stable process, Schlather’s and extremalt, and the Bayesian kriging have been put in comparison with respect to the predictive performance measures.The station by station analysis using the Bayesian approach depicted that all the 29 stations considered in this paper brought about a positive trend in the annual extreme ozone concentration. This would indicate that that the ozone concentration of all these air monitoring sites has an increasing trend with time. From the interpolation of the 5, 10, 25, 50, and 100 years’ return level, most stations found in the territories will have a higher ozone concentration.



---

*Source: 8839455-2020-11-27.xml*