Introduction

Air pollution is a complex mixture of solid particles, liquid droplets, and gases. It can come from many sources, for example: domestic fuel combustion, industrial chimneys, motor vehicles, power plants, open waste burning, agricultural activities, dust, and many other sources (Ouyang et al., 2022). According to the World Health Organization (WHO), air pollution is measured by many variables, namely PM2.5 and PM10 (particles with an aerodynamic diameter equal to or less than 2.5, also called fine, and every 10µm), ozone (O₃), nitrogen dioxide (NO₂), carbon monoxide (CO), and sulfur dioxide (SO₂) (Yang et al., 2020). Fine particles (PM10 and PM2.5) can through the lungs and then enter the body through the bloodstream, affecting all major organs (Thangavel et al., 2022; Yang et al., 2020). This can cause illness in both the cardiovascular and respiratory systems, leading to diseases such as stroke, lung cancer, and chronic obstructive pulmonary disease (Choung & Kim, 2019; Ren & Tong, 2008; Wright et al., 2023). Recent research also shows a link between prenatal exposure to high levels of air pollution and developmental delays in 3-year-old children, as well as psychological and behavioral problems later in life, including symptoms of attention deficit hyperactivity disorder (ADHD), difficulty focusing, anxiety and depression (Johnson et al., 2021; Kaur et al., 2023).

Several analytical methods are useful for monitoring air quality, including machine learning methods. Natarajan et al. (2024) applies several machine learning methods to monitor air pollution in several cities in India using k-nearest neighbor, random forest regressor, and support vector regressor models. Other analytical methods such as Markov switching models are used to understand when switching between anomalous and non-anomalous conditions in air quality occurs. The Markov switching model is a statistical analysis tool for identifying regime shifts in time series data (Franke, 2012). In the context of air quality, this can help identify periods when air pollution reaches abnormal or dangerous levels. This model allows us to group data into two or more different regimes, each with different statistical characteristics (Gao, 2020; Zakaria et al., 2019).

Numerous research studies have examined the regime-switching model specifically the Markov Switching Autoregressive (MSAR) model, applied in various field. These several studies show the versatility and effectiveness of the MSAR model in analyzing different type of data. Table 1 explaining the summary of key research employing the MSAR model, including their objective and estimation methods. These studies demonstrate the MSAR model’s ability to handle different type of data, from financial market to wind speed and mortality rates, highlighting its flexibility and robustness in various applications. The unique characteristic of MSAR is the ability to handle the regime switching without defining the threshold area first like the TAR model (Zhang et al., 2023). The original estimation method in MSAR is using EM. then the estimation method was developed using the Bayesian method to obtain optimal estimation results was proposed by Kim and Nelson (2000) and Hamilton and Raj (2002).

Table 1.

Summary Key Research Employing the MSAR Model.

There are several methods for estimating model parameters in Bayesian modeling. Some of the studies mentioned above used the Bayesian method coupled with the Gibbs sampling algorithm to estimate the parameters developed by Sims et al. (2008). Apart from that, there are also several developments in the Gibbs Sampler algorithm, including Hamiltonian Monte Carlo (HMC). HMC is an estimation method that uses the same Markov Chain Monte Carlo (MCMC) as the Gibbs sampling (Duane et al., 1987; Neal, 2011). The performance of HMC effectively mitigates the random walk behavior and correlated parameter sensitivity issues common in MCMC methods by employing first-order gradient information-based steps. The disadvantage of HMC is the number of leapfrogs, if it is too small, the algorithm will show undesirable random walk behavior, while if it is too large the computation will take a long time. To overcome that problem Hoffman and Gelman (2014) developed an extension of HMC called No U-Turn Sampler (NUTS). This algorithm can simplify the problem of HMC by automating the number of leapfrogs. NUTS also automatically stops simulation iterations if it approaches a U-Turn pattern, helping to avoid inefficient sampling of the posterior distribution.

Meanwhile, the development of the MSAR model using Bayesian estimation has been completed by Li et al. (2022) using Just Another Gibbs Sampler (JAGS) software. Their development combines the zero-inflated multilevel Poisson distribution with an autoregressive model which applied to longitudinal data. JAGS allows users to write their functions, distributions, and samplers (Wabersich & Vandekerckhove, 2014). JAGS is a development of Bayesian Inference Using Gibbs Sampling (BUGS), both of which use the MCMC algorithm for estimation. Adding new distributions to the BUGS program is complicated and requires other programs such as BlackBox Component Builder (Wetzels et al., 2010). Just like BUGS, adding new distributions to JAGS is also complicated because of testing and validation requirements, and the necessity for a clear document (Wabersich & Vandekerckhove, 2014). Different from the two previous software, in the Stan language, users can use several features to create distribution-based models that allow researchers to build based on their creativity (Annis et al., 2017). Therefore, researchers can build various models based on data-driven analysis. Modeling with Stan is widely available in several interfaces in several software. The most popular and widely used are RStan in R and PyStan in Python. With Stan available on many interfaces, it will be easier for researchers to apply the proposed method (Annis et al., 2017). In this study, we used RStan.

In real cases, not all data has a normal pattern, especially data that is suspected to have anomalous events. These anomalies can distort statistical analyses and avoid accurate modeling. In response to this challenge, researchers have developed innovative approaches such as replacing the normative Gaussian-based models with other distributions. This study approach was once carried out by Deschamps (2006) by replacing the error in the MSAR model with a Student-t distribution. Different from the approach taken by Deschamps (2006), a skewed normal Azzalini distribution which is intended by Azzalini (1985) is used to replace the error distribution carried out by Lhuissier (2019). Their studies both used the Gibbs sampling algorithm to estimate parameters. Handling skewed pattern data in this research, we propose a simulation-based model estimation using NUTS which is applied to neo-normal data distribution, namely Modified Skew Normal Burr (MSN-Burr).

This study aims to create a user-defined neo-normal Markov Switching Autoregressive Modified Skew Normal Burr (neo-normal MSAR MSN-Burr) model. Next, the proposed model is compared with the normal MSAR model to determine its ability to analyze simulation data and PM10 data. The simulation was made from three scenarios with three different distributions with the aim of finding out in general that the proposed model is able to deal with symmetric and asymmetric data conditions. The neo-normal MSAR MSN-Burr model is able to demonstrate the flexibility of the adaptive MSAR model for various data-driven distributions. The RMSE is used as the evaluation metric to better understand which model performs more effectively. Furthermore, adaptive control limits are created for each regime which are built using the highest posterior distribution for air quality mapping. A more mathematical and in-depth explanation of the MSN-Burr distribution can be seen in Iriawan (2000).

The rest of this paper is organized as follows. The next section introduces the general MSAR model. The Stan code for the general MSAR model can be seen in Osmundsen et al. (2021). The following section describes the MSN-Burr distribution and the neo-normal MSAR MSN-Burr model and demonstrates the Stan code according to the mathematical model description. The next section is comparison between the HMC and NUTS algorithm, in this section we explain the efficiency of NUTS over the HMC. After that we explain how the proposed model is estimated using a combination of EM and NUTS which we call EM-NUTS estimation. We have also provided a combination of the two estimation methods in the Stan language by adding the MSN-Burr distribution to Stan. The next section shows a simulation study comparing normal MSAR with neo-normal MSAR MSN-Burr. Simulation studies are carried out by generating data that has characteristics such as regime switching with errors using normal, double exponential, and Student-t distributions. After that we applied the neo-normal MSAR MSN-Burr model on Yogyakarta air quality data in 2021 and showing the result of the neo-normal MSAR MSN-Burr model that has converged and applied the HPD for each regime for several levels of significance. The conclusions are given in the last section.

Markov Switching Autoregressive Model

Markov switching models can be combined with time series models such as autoregressive models used to identify changes in conditions or time series data patterns (Hamilton, 1989). The forms of a Markov switching autoregressive (MSAR) model can be written in the following Equation 1.

(1)

with e_t is residual as , S_t is regime (unobserved random variable), y_t ,..., y_t-p are observation data, β₁, ... β_P are the autoregressive coefficient of order p, is the variance that is influenced by regime changes, and are mean that is influenced by changes in regime.

There are two types of modeling steps in MSAR, which are regime transition estimation and parameter estimation for each regime (Frühwirth-Schnatter, 2009). Regime transfer is an unknown condition, which is why it is called a latent variable. However, the number of regimes can be determined by various combinations. A comprehensive discussion of the MSAR model has been conducted by Kim and Nelson (1999). In this study, we refer to the normal MSAR model in Stan as a first introduction.

Neo-Normal Markov Switching Autoregressive Modified Skew Normal Burr (Neo-Normal MSAR MSN-Burr) Model

The MSN Burr distribution is a relaxation of the normal distribution developed by Iriawan (2000) from the Burr II distribution (Burr, 1942). Its cumulative distribution function (CDF) and probability density function (pdf) are

(2)

(3)

By performing the transformation, the CDF in Equation 2 and the pdf in Equation 3 are called the Modified Stable Burr or MS-Burr distribution. The mode of the MS-Burr distribution is stable at any value. However, the adjustment is needed because when compared with the standard normal distribution, N(0.1), the mode of the MS-Burr pdf value is lower. Furthermore, a transformation to fit the normal distribution such that the CDF and pdf were obtained in the form of Equations 4 and 5.

(4)

(5)

where , and σ > 0 with the µ is mode, σ is variance and lgr; is the skewness parameter of this distribution. The k value is obtained from the difference between the pdf of the standard normal distribution, N(0,1), and the pdf of the MSN-Burr distribution. Then the difference between the two distributions is equalized to zero, then the mode and scale parameter values in the MSN-Burr pdf are the same as in the standard normal pdf. The k can be written in Equation 6.

(6)

By substituting the k into Equation 5, we obtain the MSN-Burr pdf which is detailed in Equation 7.

(7)

This distribution is then added to the MSAR model in Equation 1 by replacing the residuals of the normal distribution with the MSN-Burr distribution. The definition of neo-normal MSAR MSN-Burr model is explained in Equation 8 including how to estimate its parameter.

Parameter Estimation

The combination of MSAR with MSN-Burr distribution can lead to complex estimation methods. To overcome this problem, we estimate this model using the combination of EM algorithm which proposed by Dempster et al. (1977) with NUTS which proposed by Hoffman and Gelman (2014). To enable a better understanding, we include pseudocode for both the HMC and NUTS methods, as well as a comparison of their differences. The pseudocode of HMC can be seen in Algorithm 1, meanwhile the NUTS in Algorithm 2.

Algorithm 1.

HMC algorithm.

Algorithm 2.

NUTS algorithm.

The main difference between HMC and NUTS lies in how they determine the length of the trajectory for sampling. HMC use pre-define leapfrog steps to simulate the Hamiltonian dynamics which is very good for carefully estimating parameters. In contrast, NUTS automatically define the leapfrog steps which will have an impact on trajectory length. This automatic determination is assisted by binary tree calculations which are useful for checking the U-turns. The addition of the automatic leapfrog determination step and the binary tree calculation in the NUTS algorithm prevents it from getting trapped in local optima, enabling it to achieve optimal parameter estimates more efficiently.

Combination EM-NUTS Estimation

The key combination of the EM and NUTS is only using the Expectation step in EM then the Maximization step using NUTS. Before we calculate the expectation step, we need to find the likelihood using the joint density of y_t, s_t and s_t-1.

1. Deriving the joint density function of y_t, s_t and s_t-1 conditionally on past information ψ_t-1

Where Ψ_t-1 is the observation value on t − 1 of time, while yt is the observation value at the time t which follow the MSAR model using MSN-Burr distribution on Equation 7 with parameter is given by Equation 8.

(8)

Since the s_t parameter was added to the MSN-Burr distribution, the MSN-Burr MSAR model equation has parameters in each regime. Mode is µ_st, variance is , skewness is , autoregressive parameter is .

2. Find the function ) by integrating s_t and s_{t − 1} from the joint density over all possible values of s_t and s_{t − 1}

where is the marginal density which is the average of conditional density weighted by

where M is the number of regimes inside the model. The likelihood function is written below.

To calculate is going to be solved using the filtering process. In short, the filtering process aims to get the value.

After the Expectation step, we move to Maximization step using the NUTS algorithm. In this step, we aim to optimize the parameter .

1. Initialize the parameter value in each regime θ _st and step size ε on the first iteration.

2. Draw initial momentum of at every iteration from

3. Use the leapfrog integration method to simulate Hamiltonian dynamics. There are three main steps:

   • a. Update momentum (half-step):

     

   • b. Update position:

     

   • c. Update momentum (half-step):

     

4. Build the binary tree

   • a. Expand both forward and backward in time:

   • - Forward:

     

   • - Backward:

     

   • b. Check the U-Turn

     

5. Calculate the acceptance rejection

   • a. Calculate the Hamiltonian for the initial and proposed states:

     

   • b. Calculate the acceptance probability

     

   • c. Accept or reject the proposed θ_st:

     

6. Repeat the procedures for the desired number of iterations to produce samples from the specified distribution

The NUTS algorithm is computationally intensive because it requires calculating the binary tree structure at each iteration in order to dynamically adjust the trajectory length. This process involves multiple leapfrog integration phases for both forward and backward expansion, as well as repeated checks for U-turn condition. This algorithm will make the parameter estimation process in the MSN-Burr MSAR model in Equation 8 converge more quickly. We briefly summarize the overall combination of the two estimation methods which can be seen in Algorithm 3.

Algorithm 3.

Combined EM-NUTS estimation method.

To simplify the estimation process using the proposed algorithm, we use the Stan language available in R. The Stan language was developed to overcome convergence problems that commonly occur in Bayesian inference using Gibbs sampling (BUGS) languages (Gelman et al., 2015). The key stages of the Stan Language are data and model input, calculating the log of the pdf and their gradients, a warm-up phase for parameter tuning, applying NUTS, monitoring convergence, and calculating summary inference (Hoffman & Gelman, 2014). Standard distributions such as normal, Poisson, binomial, etc. are already available in the booth. However, Stan is very flexible in adding new distributions by writing a log of pdf from that distribution. Stan uses a numerical auto-differentiation method by utilizing reverse-mode automatic differentiation to automatically perform function reduction (Carpenter et al., 2017).

Simulation Study

From Equation 7, when the value of the skewness parameter λ = 1, it can be seen that the MSN-Bur distribution becomes a normal distribution, which is symmetric and bell-shaped. Therefore, normally distributed data are used to validate the user-defined MSN-Burr distribution in the Stan program, allowing Stan to use the user-defined MSN-Burr distribution for estimation. The fact shows that MSN-Burr can accurately and consistently estimate µ and σ in accordance with the generator parameters serves as evidence of its validity in detecting normal data.

We then present comparative evidence between the standard conventional MSAR model and the neo-normal MSAR MSN-Burr model in Stan. We created three different scenarios involving the normal distribution, the double exponential, and the Student-t distribution. Three distributions are used to demonstrate that neo-normal MSAR MSN-Burr is capable of collecting symmetric distribution characteristic data. These scenarios contain a challenging scheme, which is used to compare the capabilities of the neo-normal MSAR MSN-Burr model with the standard MSAR normal model to detect data with different variances and zero-centered data using leptokurtic properties. In each scenario, we generate two different time series observations based on the selected distribution. As many as 700 observations are generated based on the selected distribution, then 300 observations with different parameters to obtain a total of 1,000 observations with a proportion of 0.7 in regime 1 and 0.3 in regime 2. We then estimate each regime in each scenario using a standard autoregressive model to find the target parameter. A visualization of the generated non-normal time series data is displayed in Figure 1. The target parameters are then derived based on the autoregressive models of each scenario and are available in Table 2. Target parameter in Table 2, we assign a value to the slope parameter λ = 1 to mimic the symmetric conditions that would be estimated by neo-normal MSAR MSN-Burr.

Figure 1.

Time series plot of generated data using (a) normal, (b) double exponential, and (c) Student-t distributions.

Table 2.

The Target Parameter from Three Scenarios for Normal, Double Exponential, and Student-t Distribution.

From the three scenarios above, we obtain the 95% credible interval for each parameter as seen in Table 3. We can use the 95% credible interval in Table 3 as a parameter significance test by providing a range of values for the parameter of interest that are considered plausible given the observed data. In the context of hypothesis testing, if the target parameter value in Table 3 falls within the interval, it suggests that the data is consistent with that value, supporting the null hypothesis. Conversely, if the target parameter is outside the credible interval, it leads to reject the null hypothesis (Martinez & Martinez, 2016).

Table 3.

The 95% Credible Interval for Estimated Parameters of Three Scenario Simulations for Normal, Double Exponential, and Student-t Distributions.

Each scenario was replicated 200 times. In each replication, estimation is carried out using the NUTS algorithm which is available on Stan through 10,000 iterations. Since each scenario was replicated 200 times and in every one of them was calculated using the sampler from NUTS, we computed the root-mean-square error (RMSE) to evaluate and compare the quality of parameter estimates within each scenario. The RMSE for estimated is defined in Equation 9 (Walther & Moore, 2005):

(9)

where θ is the list of target parameter from Table 2, is the estimation result in each replication, and n_rep is the number of replications. Since the θ is the list of target parameter, the RMSE will be calculated for each parameter. RMSE performance can be seen from a value close to zero. When the RMSE approaches zero, the parameter estimate is better.

All the target parameter in Table 3 lies within the 95% credible interval. Likewise, for that the MSN-Burr distribution in the neo-normal MSAR MSN-Burr model is able to estimate data generated from symmetric distributions, namely normal, double exponential, and Student-t distributions as can be seen at the lgr; value in each regime that is close to one. For the autoregressive parameters, namely β₁₁ and β₁₂, they have their own hypothesis, namely that they must be inside unity, which means they do not have a zero value in the credible interval. If the autoregressive parameter is zero then the model essentially becomes a random walk. The evidence that none of the intervals for the parameters β₁₁ and β₁₂ have a value of zero is presented in Table 3. So, it can be concluded from the results of the 95% credible interval the neo-normal MSAR MSN-Burr model has succeeded in identifying random data points with a symmetrical distribution.

The results of the RMSE from each parameter calculation are shown in Table 4. This RMSE was calculated through 200 replications, with each replication generate estimated value from each parameter. This approach allows us to assess the error for each parameter across all replications. In the generated normal distribution scenario, the RMSE value of the sigma parameter in the neo-normal MSAR MSN-Burr model is slightly higher than in the normal MSAR model. In the generated double-exponential and generated Student-t distribution scenarios, the RMSE value of the p parameter in the neo-normal MSAR MSN-Burr model is slightly higher than normal MSAR. This doesn’t mean that the neo-normal MSAR MSN-Burr has poor performance. From several explanations given, it can be concluded that the MSAR model that uses MSN-Burr distributed error can detect symmetric random data points from the normal, double exponential, and Student-t distribution.

Table 4.

Root-Mean-Square Error (RMSE) for the Estimated Parameter of Three Scenario Simulations for Normal, Double Exponential, and Student-t Distributions.

We also calculate the RMSE for each observation to find out whether the predicted value is around what the observation should be or not. This RMSE comparison was carried out by the normal MSAR model with the MSN-Burr MSAR to find out which model is better at capturing symmetrical patterns in the three simulation scenarios. RMSE calculations for predicted values from observations are presented in Table 5.

Table 5.

Root-Mean-Square Error (RMSE) for Predicted Observation of Three Scenario Simulation for Normal, Double Exponential, and Student-t Distributions.

The observation RMSE calculation produces a percentage value of the difference between the normal MSAR model and the neo-normal MSAR MSN-Burr in the three scenarios. Table 5 shows that neo-normal MSAR MSN-Burr is 14.0408% better in modeling the generated normal data scenario and 82.7683% in modeling the generated student-t data compared with the normal MSAR model. However, both MSAR models seem unable to model data that has very high leptokurtic properties in the Generated Double Exponential scenario. The percentage difference in this scenario is only 0.0046%, which shows that the two models provide the same estimation results for estimating the generated data with high leptokurtic properties.

Application

This session discusses the application of the neo-normal MSAR MSN-Burr model using an air quality dataset from Yogyakarta, Indonesia, collected in 2021 and recorded hourly. This dataset is openly available from the Yogyakarta City Environmental Service. This dataset has 8,760 observations and has several missing values, we handle it by employ moving average. We use this method by taking the average of the previous 24 observations, because the data used is data taken hourly so the assumption of taking the average of the previous 24 hr.

As in the introduction, there are several variables in air quality, namely PM2.5 and PM10 (particles with an aerodynamic diameter equal to or less than 2.5, also called fine, and every 10 µm), ozone (O₃), nitrogen dioxide (NO₂), carbon monoxide (CO) and sulfur dioxide (SO₂) (Manisalidis et al., 2020). In this study, we use fine particle PM10 variable to demonstrate the neo-normal MSAR MSN-Burr model. The detailed PM10 movement patterns are recorded hourly and its histogram which looks skewed to the right can be seen in Figure 2.

Figure 2.

(a) Time series plot of PM10 and (b) histogram of PM10.

This application compares the standard PM10 level published by WHO with the actual accidents that occurred in Yogyakarta city in 2021. Clean air typically contains PM10 levels of less than 15 micrograms per cubic meter (µg/m³) (World Health Organization, 2021). The results obtained from the neo-normal MSAR MSN-Burr model are in two states: normal and abnormal.

Results and Discussions

This section shows the results obtained from the neo-normal MSAR MSN-Burr model written in the Stan code available on Github  https://github.com/DwilaksanaAbdullahRasyid/MSAR_MSN-Burr. The normal distribution was chosen as the prior over several parameters because it can approximate many other distributions. This makes it a good choice for modeling a wide range of phenomena (Kruschke, 2015). Bayesian methods are very effective when used on small amounts of data. Indeed, Bayesian uses Monte Carlo simulation (Beer et al., 2013). Therefore, if this estimation method is applied to big data, including air quality data, a computer device with large capacity random access memory (RAM) is required. This happened in this study when the neo-normal MSAR MSN-Burr model parameter estimation for the PM10 variable on the Yogyakarta 2021 air quality dataset was implemented on a computer with a RAM capacity of only 12 GB and the sampling scenario was repeated four times. The program is run for 100,000 iterations in each replication, consuming 11 GB of memory, so it is difficult for the parameters to reach convergent conditions. To overcome this problem, a sampling scenario is applied by reducing the domain range of each parameter only in the predicted area with a feasible pdf in the next replication. This scenario provides excellent parameter estimation results, it is found that the slope parameter of the MSN-Burr distribution shows a dominant role in achieving convergence. The membership of each observation after all parameters converge is shown in Figure 3, with regime 1 representing PM10 under normal conditions and regime 2 representing PM10 under abnormal conditions.

Figure 3.

Membership distribution from observations (a) for regime 1 and (b) for regime 2.

In the fourth repetition, the parameters converge, as indicated by the perfect “Rhat” value of all estimated parameters equal to one. “Rhat” as a tool to monitor MCMC convergence was discovered by Gelman and Rubin (1992). This value measures the ratio of the average variance of the samples in each chain to the variance of the pooled samples across the chain. According to Gelman and Rubin, independent Markov chains should be sampled until all values of “Rhat” are equal to one after being initialized with diffuse initial values for the parameters. According to Gelman and Rubin, independent Markov chains should be sampled until all values of “Rhat” are equal to one after being initialized with diffuse initial values for the parameters. Some of these parameters are for the probability regime shifts, for the intercepts, for the autoregressive constants, for the variances, and for the skewness. A summary of the model estimation results can be seen in Table 6 and the estimation results of the MSAR normal model using original EM algorithm as a baseline can be seen in Table 7.

Table 6.

The Estimated Parameter of Neo-Normal MSAR MSN-Burr with Two Regimes.

Table 7.

MSAR Normal Estimated Parameter.

The MSAR parameters presented in Table 7 indicate that the estimate for the autoregressive parameter, β₁₂, is 1. This estimation suggests that β₁₂ parameter leads to a random walk, which implies that the model’s predictions do not revert to a mean but instead exhibit continuous and irregular drift over time. As a result, regimes are separated irregularly across the model, indicating a high level of uncertainty and unpredictability in regime changes.

To get a deeper understanding of the impact of the random walk on the model, we may look at the regime shifts in the conventional MSAR model, as seen in Figure 4. The MSAR normal regime changes are contrasted to those of the neo-normal MSAR MSN-Burr model. The main difference is the separation of observations within each regime. The separation of the MSAR normal model in Figure 4(a) is highly irregular, showing unpredictability between regimes. In contrast, the proposed model in Figure 4(b) shows the separation of regime and can explain skewness in the data. The separation to regime 1 explaining in-control conditions, whereas regime 2 explaining out of control conditions. This separation demonstrates the capability of the neo-normal MSAR MSN-Burr in distinguishing different regime, leads to reliable and interpretable results.

Figure 4.

Regime membership (a) MSAR normal (b) neo-normal MSAR MSN-Burr.

The performance of the models is compared using the Root Mean Square Error (RMSE) metric. The normal MSAR model has an RMSE of 11.87909, whereas the MSN-Burr MSAR model achieves a significantly lower RMSE of 9.9303. With the resulting RMSE the normal MSAR is 11.87909 and the MSN-Burr MSAR is 9.9303. This significant reduction in RMSE with 16.4049% improvement illustrates the proposed greater accuracy and effectiveness in capturing underlying data patterns when compared to the MSAR normal model.

After estimating using EM-NUTS, calculate the Upper Control Limit (UCL) and Lower Control Limit (LCL) using the Highest Posterior Distribution (HPD). HPD is a method for estimating confidence intervals from asymmetric distributions (Gelman et al., 2014). Chen and Shao (1999) estimated HPD intervals using Monte Carlo procedures. The basic idea for creating a (1 - α) HPD interval (also called as a credible interval) is to use the concept of density equilibrium by solving two simultaneous equations. For example, suppose that f (x) is a pdf of random variable X. The values a and b are the lower and upper limits of the HPD interval in the domain so that . Then the values of a and b respectively as the lower limit and upper limit of the credible interval must cover a certain area of the pdf, so that , for a certain level of significance α. This concept can be stated that the choice of a and b which maintains both the highest density and the area under the density between the chosen boundaries, can be obtained by solving the following two equations simultaneously.

(10)

(11)

In general, the calculation of ULC and LCL on HPD can be seen in Algorithm 4.

Algorithm 4.

Control Limit Estimation Algorithm in HPD.

The HPD calculation is applied to the MSN-Burr distribution in each regime. The aim adding the HPD in this study is to know the early warning limits before a regime change occurs. The results of the control limit estimation using Algorithm 2 are the upper and lower limits of each regime. Several significance levels were tested to be applied to Algorithm 2. We use three different levels of significance to see how sensitive the specified control limits are. The significance levels used are 1%, 5%, and 10%, as shown in Table 8.

Table 8.

Several Levels of Significance in Regime 1 and Regime 2.

Figure 5 provides a visualization of the boundaries of each regime using HPD with three different significance levels. As described at the beginning of Application section, regime 1 is a normal PM10 condition in Figure 5(a), and regime 2 is an abnormal PM10 condition in Figure 5(b). In addition to the visualization, Figure 5 also shows the different modes for each regime. The normal condition mode of PM10 in Yogyakarta is 8 µg/m³ and the abnormal condition mode is 19 µg/m³.

Figure 5.

HPD for several significance level, red line for 10%, green line for 5%, and blue line for 1%, which are overlaid with a histogram for each regime (a) regime 1 and (b) regime 2.

High levels of air pollution, particularly the PM10 variable in the air quality measurements, can be seen in Figure 5(b), especially under abnormal conditions. High levels of PM10 can cause health problems, as described in Introduction section. Therefore, a large level of significance is required to create highly sensitive control limits. Under normal conditions, however, all of the upper control limit compliance the WHO air quality guideline. This approach aligns better with WHO guidelines and ensures more sensitive detection of deviations. To address abnormalities, adaptive measures such as real-time monitoring adjustments are essential for maintaining air quality within safe limits.

Conclusions

The MSAR model was created to capture the movement between certain conditions and other conditions. However, with anomalies in real-world conditions that occur successively, the standard MSAR model cannot capture them because the normal residual assumption is violated. Therefore, the proposed model contains residuals with a neo-normal distribution, namely MSN-Burr. With the help of the Stan programming language using EM-NUTS as an estimator, we obtain a new model that is mathematically complex but capable of being developed.

This paper presents the simulation study using normal, double exponential, and Student-t distribution. All three have something in common, which is a symmetrical pattern. The 95% credible interval results, the neo-normal MSAR MSN-Burr model can estimate data generated from the normal, double exponential, and Student-t distributions which can be seen lgr; value in each regime that is close to one. We then also provide bias and RMSE calculations in each parameter to compare the neo-normal MSAR MSN-Burr model and normal MSAR. As a result, the neo-normal MSAR MSN-Burr model can capture the conditions of symmetric pattern data generation.

This paper also explained the calculation of RMSE from predicted value in simulation study. The double exponential scenario, the RMSE improvement percentage of neo-normal MSAR MSN-Burr is not much different from MSAR normal. This is because the double exponential distribution generates a high leptokurtic pattern. However, in the other two scenarios, namely the normal and Student-t scenarios, neo-normal MSAR MSN-Burr has better performance than MSAR normal in producing RMSE.

The neo-normal MSAR MSN-Burr model can capture shifts in PM10 conditions, where each condition exhibits an asymmetric pattern. PM10 is a type of particle with an aerodynamic diameter of 10 µm, which can penetrate the lungs and then enter the body through the bloodstream, affecting all major organs. This can cause diseases of both the cardiovascular and respiratory systems. In Yogyakarta, the PM10 in normal conditions is 8 µg/m³ which is classified as slightly higher than WHO standards.

In addition, we also calculate the upper and lower limits for each condition using HPD. All of the significance level are advised to be implemented in Yogyakarta due to its lower limits more settle to the WHO air quality guideline. The higher significance level is produces control limits that are more sensitive than the lowest significance level. This approach could detect the abnormalities and the significance level which compliance to the guidelines.