Infrastructure deterioration modeling with an inhomogeneous continuous time Markov chain: A latent state approach with analytic transition probabilities

Markov chains have been widely used to characterize performance deterioration of infrastructure assets, to model maintenance effectiveness, and to find the optimal intervention strategies. For long‐lived assets such as bridges, the time‐homogeneity assumptions of Markov chains should be carefully checked. For this purpose, this research proposes a regime‐switching continuous‐time Markov chain of which the state transition probabilities depend on another, latent, Markov chain that characterizes the overall aging regime of an asset. With the aid of a state‐augmentation technique, closed‐form solutions for the transition probabilities are analytically derived, making the statistical analysis simple. A case study is presented using the open Ontario Bridge Condition data for provincial highway bridges. The case study demonstrates that the proposed method allows to (1) estimate a statistically superior model to the homogeneous Markov chain and (2) obtain results with comparable accuracy in approximately 48% of the computation time of the state‐of‐the‐art inhomogeneous Markov chain.


INTRODUCTION
Deterioration modeling is a key element of an infrastructure management system Prakash et al., 2020). Markov chains have been widely used to model infrastructure deterioration processes as a stochastic transition process among deterioration states. The popularity can be attributed to at least the following three aspects of the Markov chain model: First, for a discrete-time homogeneous Markov chain, which is the most popular one among all, both parameter estimation and model prediction are simple and straightforward. The mathematical tractability is also a very attractive aspect in practice. Second, maintenance effectiveness is another important aspect of deterioration modeling (Chu & Durango-Cohen, 2008;Zhang & Gao, 2011) that can be readily characterized and modeled along with the natural deterioration within one single Markov chain model (Lin et al., 2019). Third, with a Markov chain deterioration model, the optimal management policies can be formulated as a Markov decision process and solved using a dynamic programming method. Using different sets of decision variables, the Markov decision process can be used to conveniently estimate required future interventions and the corresponding expenditures both at an aggregate level (e.g., the early pavement management system developed by Golabi et al., 1982) and at a segment level (most recent applications belong to this category). Therefore, although many advanced deterioration prediction models have been developed (e.g., artificial intelligence), Markov chains are still often used to characterize deterioration processes in order to find the optimal intervention policies for various infrastructure systems, e.g., Gao and Zhang (2013) and Mizutani et al. (2020).
Over the years, a number of enhancements have been introduced to relax the assumptions embedded in the Markov chain deterioration models. These include, for example, the extensions from discrete time to continuous time, from homogeneous to inhomogeneous (or nonstationary) transitions, and from the Markov chain to the semi-Markov chain. In addition, several computational challenges related to statistical estimation and deterioration prediction have also been addressed. More details will be reviewed in Section 2. Despite these, the modeling of time-inhomogeneity and long-range dependence exhibited in the performance deterioration of infrastructure assets remains to be an important challenge that has not been adequately researched. Although the traditional inhomogeneous Markov chain can be used to model the apparent time-varying transition rates from one state to another state, there are other time-dependent phenomena that invite more sophisticated modeling techniques. For example, Robelin and Madanat (2007) describe a threedimensional Markov chain model used to characterize long-range history dependence that the traditional onestep Markov chain would not be able to appropriately model. From another perspective, Bocchini et al. (2013) explain a situation where the transition probabilities from a serviceable state to a maintainable state for a new bridge are different than those of a bridge after maintenance and/or rehabilitation. As a result, for two infrastructure assets-one immediately after construction and the other after major maintenance or rehabilitation-although both exhibit the same "intact" state, their probabilities of moving to another advanced deterioration state are different. If the timing and the type of specific maintenance actions were properly documented, the maintenance effects on the transition probabilities or rates could be explicitly modeled by casting them as a function of the maintenance activities. In practice, however, the historical maintenance records are often unavailable, rendering this modeling approach very hard to implement, if not completely untenable (Lin et al., 2019). In addition, when the working environment defined by those factors change from one set to another, the whole deterioration trend may change. Unfortunately, and similar to the case of maintenance history, the changing working environment is often not carefully characterized and recorded either. As a result, many time-varying factors or physical mechanisms that affect asset performance deterioration rates are either not available or unobservable. Therefore, it is hard to characterize the effects of the unob-servable factors on the deterioration using the traditional, regression-type modeling approach.
For this reason, it is desirable to use a model that allows for modeling the hidden changes in the transition rate over time to assess whether there is a change in the transition rate over time and, if so, to forecast the future deterioration process by taking this change into account. Meanwhile, it is also important to specify the Markov chain model so that the process of transition rate change over time is not overly complicated, computational efficiency is ensured, and appropriate information is provided as input for decision making (e.g., Markov decision process).
In this study, an unobservable potential deterioration stage is defined and modeled as a variable referred to as latent state (LS), which stochastically changes over time in accordance with a homogeneous discrete-state continuous-time Markov chain (CTMC), and conditional on a given LS, the progression of the observable deterioration condition (DC) is characterized by another homogeneous CTMC. The resultant deterioration process, therefore, exhibits an inhomogeneous property. Specifically, the transition rates of DCs are different even for the same DC if the LS is different, and the LS can be interpreted as an indicator of the stage of deterioration (e.g., early, middle, or end stage) of the infrastructure asset. The proposed model can be interpreted as a hidden Markov model or a regime-switching Markov chain. As it will be shown later in the paper, the proposed model is advantageous to other inhomogeneous models, in that, its transition probability has an analytic solution. This is the most important originality of this research, and to the best of the authors' knowledge, there have been no previous studies that utilize hidden Markov models or regime-switching models (i.e., define the LS) to derive analytic solutions for transition probabilities. In addition, due to the Markovian property, the proposed model can be readily integrated into a Markov decision process or a partially observable Markov decision process for maintenance policy optimizations.
Note that since the objective of this study is to develop a methodology that contributes to asset management processes of budget planning and repair quantity estimation from a long-term perspective (e.g., processes in "6 Planning" of ISO 55001 International Organization for Standardization, 2014), a method is proposed to estimate the deterioration process once using all the inspection data currently held by the asset owner. In such processes, updating of the deterioration process by additional observation of new inspection data is not covered, because excessively repeated updating of the deterioration process is considered undesirable for long-term planning. Of course, updating the deterioration process in real time or at shorter time intervals is required in other processes of the asset F I G U R E 1 Overall framework management system (e.g., operational processes), but is beyond the scope of this study.
The remainder of this paper is organized as follows. Section 2 includes a literature review. In Section 3, the proposed model is set up, and the closed-form solution of the transition probabilities is derived. Section 4 introduces a statistical estimation method for the proposed model by formulating the likelihood function. After that, a simulation example is presented in Section 5 to confirm the identifiability of the proposed model, followed by a case study using a real-life bridge condition data set in Section 6. Finally, Section 8 concludes the study with major findings and suggestions for future work. The overall framework of this study is shown in Figure 1.

LITERATURE REVIEW
There is a huge body of literature related to Markov models and their applications to infrastructure deterioration modeling. It is beyond our capacity, neither is it our intention, to provide a comprehensive review here. Rather, we try to provide a high-level brush over the historical development of Markov deterioration modeling from the assumption relaxation and numerical challenges perspective. The whole point is to show readers that the latent inhomogeneity in the transitions between DCs was only a recent research topic after all the apparent deterioration phenomena have been well researched, and thus that the proposed LS modeling approach itself is truly original and novel.
Among all stochastic deterioration processes (excluding traditional regression models), Markov chain was the first model that was applied in actual infrastructure management systems. As a matter of fact, the very first generation of computerized pavement and bridge management systems both used the Markov chain (Golabi & Shepard, 1997;Golabi et al., 1982). The first monograph that systematically introduced the Markov chain and its application to deterioration (mainly metal fatigue) was due to Bogdanoff and Kozin (1985). The model was later seen to be applied to other infrastructure systems such as sewer pipes (Abraham & Wirahadikusumah, 1999;Baik et al., 2006), stormwater pipes (Micevski et al., 2002), nuclear piping (Fleming, 2004), building facilities (Jin & Mukherjee, 2014), and offshore structures (Zhang et al., 2017). The following review will mainly focus on the applications to pavement and bridges.
Over the past four decades, the Markov chain model has received a lot of attention, and several enhancements have been introduced to improve its flexibility and prediction accuracy. First of all, the homogeneity assumption that is commonly adopted in the Markov model was reviewed and extended. For pavements, Butt et al. (1987) employed a staged-homogeneous Markov model in which the transition probabilities were defined for each stage that changes deterministically over time, whereas Li et al. (1997) used an inhomogeneous discrete-time Markov chain to predict asphalt pavement performance, with transition probabilities determined by Monte Carlo simulation. For bridge deck deterioration, Morcous (2006) was probably the first researcher in infrastructure management who for the first time seriously examined the effects of time inhomogeneity on the accuracy of deterioration prediction.
The aforementioned discrete-time Markov chain models and many others alike (Chiachío et al., 2020;Kosgodagan-Dalla Torre et al., 2017) require the deterioration data be time series. To adapt to nonperiodic inspections or missing deterioration data often experienced in practice, CTMCs (Ferreira & Suchard, 2008;Geweke et al., 1986) were used. An extra benefit of the CTMC is that the inhomogeneity in the deterioration process can be flexibly modeled through the use of transition rate functions, rather than constant rates (Gokhale et al., 2004). For example, Kobayashi et al. (2010) proposed to use the Weibull hazard function to represent the time inhomogeneous CTMC and derived its transition probabilities in nonanalytic solution form.
A challenge related to the inhomogeneous CTMC is the solution of the Kolmogorov-Chapman equation, which now is written in its differential form. The normalization technique is often used to numerically analyze time inhomogeneous CTMCs (Arns et al., 2010). For general time inhomogeneous CTMCs, where the dynamic change process of the transition rate is not specified, Li et al. (2013) describes four types of numerical solution techniques: Runge-Kutta method, uniformization, Monte Carlo simulation, and state-space enrichment.
The next assumption examined was the conditional independence assumption, for which two extensions were introduced. One was the idea of a semi-Markov chain in which the probability of transition from one state to another depended not only on the pair of states before and after the transition but also on the duration of the preceding state being stayed in (Ariza et al., 2020;Black et al., 2005;Manafpour et al., 2018;Nesbit et al., 1993;Thomas & Sobanjo, 2016;Yang et al., 2009). The other approach taken to relax the one-step conditional independence assumption was due to Robelin and Madanat (2007), who proposed a history-dependent Markov model using a state-augmentation technique. Several comparative studies have been conducted to investigate the pros and cons of homogeneous and inhomogeneous Markov chains as well as the semi-Markov chains (Ariza et al., 2020;Yamany et al., 2021). Related to the semi-Markov model, Wellalage (2015) proposes a methodology to posteriorly estimate a Weibull hazard model from the estimation results of a homogeneous Markov chain.
The further next enhancement deals with measurement error and missing data, a practical issue in empirical deterioration model. Kobayashi et al. (2012) and Lethanh and Adey (2013) proposed "hidden Markov chain" models, in which the observed condition states were considered to be a realization from the true (but directly unobservable) Markov chain contaminated by "specification bias." To deal with missing records before and after maintenance actions for the modeling of maintenance effectiveness, Lin et al. (2019) proposed an advanced statistical estimation method that integrates a data augmentation technique into the Markov chain Monte Carlo (MCMC) simulation. Calvert et al. (2020) proposed a simple maximum likelihood estimation procedure to deal with truncated inspection data for railway bridges. In addition, the past 5 years have also witnessed a few interesting advancements, including Bayesian model updating (Mishalani et al., 2018), perturbed Markov process (Oumouni & Schoefs, 2021), and statistical estimation using convolutional neural networks (Liu et al., 2022).
Despite the large number of previous studies, to the best of the authors' knowledge, no models have been developed to characterize the latent inhomogeneity described in Section 1. As discussed above, the "hidden Markov model" developed by Kobayashi et al. (2012) was intended to account for measurement errors. Lethanh et al. (2015) proposed a "Poisson hidden Markov model" but it was to capture the interaction between the occurrences of potholes (modeled as a Poisson process) and surface roughness (modeled as a Markov chain). In addition, both studies did not seek an analytic solution of the transition probability; rather, they used the completed likelihood function to simultaneously sample both the parameters and the random latent variables. To this end, the continuous time axis had to be discretized into small time intervals, in which the LSs were then sampled, adding a lot of computational loads to the statistical analysis.
Therefore, a proper approach to modeling the hidden inhomogeneity in the deterioration rate is still lacking. However, in other research fields, particularly in the structural health monitoring field, a regime-switching modeling approach has emerged to be a recent development for the modeling of hidden or latent changes, for example, Zhang et al. (2021) and Wu and Ding (2022). Still, none of these explored the possibility of a Markov chain whose transition rates were modulated by another CTMC. This is the type to be explored in this study.

State and transition rate
A deterioration process of an infrastructure asset is modeled. The notation used hereafter is listed in Appendix A.
The concept of the proposed model is illustrated in Figure 2, which can be referenced along with the following formulation. The DC state space is defined as  = {1, 2, … , }, for which the larger the value of DC , the more deteriorated the asset. Variable ∈  = {1, 2, … , } is defined to represent the LS of the transition rate. The larger the LS, the greater the transition rate, indicating a later stage of deterioration. Assuming no measurement error and no spontaneous recovery of the DC, the transition rate matrix of the DC can be defined as an upper triangular matrix, which represents the probability that the DC transitions from 1 ∈  to 2 < 1 is 0. This assumption is commonly set when modeling deterioration processes by Markov chains (Madanat et al., 1995). This assumption could be relaxed (Kobayashi et al., 2012), but it is not considered here to prevent the model from becoming overly cumbersome. Therefore, when applying the proposed model to inspection data, which may contain excessive measurement error, attention should be paid to estimation bias. The transition rate matrix (an infinitesimal rate matrix) , which depends on the LS , is set as , , ∀ ∀ , and , , = 0 as is the absorbing state. To avoid excessive model complexity and ensure the identifiability of the parameters, it is assumed that the LS transition does not depend on the DC and that the transition probability of the LS from 1 ∈  to 2 (< 1 ) is 0. Although these assumptions might theoretically be relaxed, they are not considered in this study. Based on these assumptions, the transition rate matrix of the LS is set as an upper triangular matrix independent of the DC as where , = − ∑ = +1 , ∀ , and , = 0. In and , , 1 , 2 ( ∈ ; 1 ∈ ; 2 ∈ { 1 + 1, … , }) and 1 , 2 ( 1 ∈ ; 2 ∈ { 1 + 1, … , }) are nonnegative parameters that represent the deterioration rates of the DC and LS, respectively, and are statistically estimated from the observed data. The larger the value of these parameters, the faster the deterioration.
The state of the asset, ∈  = {1, 2, … , }, is represented by a combination of the DC and the LS and is defined as a state-augmentation operation: The novelty and advantage of this study compared to the other existing studies are that it derives a timeinhomogeneous CTMC whose transition probability has the analytic solution by (1) defining the LS in addition to the DC and expressing the time transition of the LS as a CTMC and (2) defining the state by the state-argumentation operation in Equation (3). The use of analytic solutions for the transition probabilities allows the likelihood or posterior probability density to be expressed analytically, which reduces the computational load by avoiding numerical integration as shown in Section 6.3.2, and also improves compatibility with gradient-based algorithms such as the Newton-Raphson method and the Langevin Monte Carlo. Another advantage of the proposed model is that the transition of the state is represented by a time-homogeneous CTMC, which has the universal characteristics of the traditional CTMC, as described in Section 3.2. The transition rate matrix, which represents the time evolution of the state, can be represented as a block matrix as and is the × identity matrix. By defining Equations (1), (2), and (3), the transition rate matrix of the states is uniquely determined, and in the case of this study, and also become upper triangular matrices. This means that the state 1 ∈  does not transition to 2 (< 1 ).
Here, let us explicitly state the meanings of the LS, DC, and State. In this study, these three terminologies are defined as follows: LS: The LS is the state variable that represents the unobservable internal deterioration mechanisms that affect the deterioration rate. The essence of the definition of the LS is that it is unobservable.
Although the existence of the LS is a prerequisite for the proposed model, it can also be ascertained and interpreted a posteriori from the statistical estimation results of the model. In the general time-inhomogeneous models, the transition rate is represented as a continuous function, but the transition probability does not have an analytic solution.
On the other hand, the proposed model defines the LS and represents its transition process as an exponential hazard model, so that the transition probability has an analytic solution, thereby representing the time-inhomogeneous deterioration process in an alternative and simplified manner. DC: The DC is a state variable that represents the revealed or observable state of infrastructure deterioration as defined in the inspection guidelines. The DC is the same for the same state of infrastructure deterioration, regardless of the value of the LS. The essence of the definition of the DC is that it is observable. State: The state is an infrastructure deterioration state variable defined by the augmentation of the LS and the DC. Since the LS is unobservable and the DC is observable, the state is partially observable. For example, if = 3 and = 2, observing DC 2 means that the state is 2 (= 2 + (1 − 1) × 3) or 5 (= 2 + (2 − 1) × 3). This property, which restricts the possible states based on the observed DC, is referred to here as partial observability. The state is defined to represent the time-inhomogeneous transition of the DC as a time-homogeneous CTMC of the state, and mainly plays a role for the mathematical formulation.
The LS transition itself is unobservable, but can be indirectly observed as a difference in the transition time between the same DCs (e.g., if the observed time that the DC transitions from 1 to 2 is long, the LS is estimated as 1; if it is short, the LS is estimated as 2). From observed data containing only the time transition of DCs, both (1) the transition rate matrix of the DC, and (2) the transition rate matrix of the LS, , can be estimated.

Markov transition probability
The Markov transition probability matrix of the state evolution, ( ), is defined as where > 0 is the time interval. ( ) depends only on the current state , that is, the combination of DC and LS , not on the past transition processes. In practice, LS is unobservable. Therefore, for practical use, in a situation where the current DC and the elapsed years since construction are observed, the following procedure can be used to analytically derive the occurrence probability of DC at a future time point: (1) Calculate the occurrence probability of the current LS based on the elapsed years, (2) Calculate the occurrence probability of the future state based on the occurrence probability of the current LS , and (3) calculate the occurrence probability of future DC by marginalizing the occurrence probability of future state with respect to LS . In addition, ( ) satisfies the time consistency condition with respect to any positive real number as The transition probability can be found by solving the Kolmogorov forward equation (a first-order differential equation): It is known that the solution of Equation (10) is given by a matrix exponential (Ferreira & Suchard, 2008): The eigenvalues of can be found as ( 1,1 , 2,2 , … , , ) = ( 1,1 , 2,2 , … , −1, −1 , 0), where , ( , ∈ ) is the × th element of . Assuming that , ≠ , ∀ , ∈ , the characteristic polynomial of has distinct roots, so is diagonalizable and can be factorized as where = diag( 1,1 , … , , ), and is the × invertible matrix whose th column is the eigenvector , of , which corresponds to the eigenvalue , . The matrix exponential in Equation (11) can then be easily solved, and the Markov transition probability matrix for an arbitrary time interval (> 0) can be derived as where e, = diag(exp( 1,1 ), … , exp( , )). Based on Equation (13), the Markov transition probability , ( ) can be derived as where , , and , , are × elements of and −1 , respectively. Note that in Equation (14), the following two conditions are automatically satisfied: , ( ) = 0 if > and ∑ = , ( ) = 1. In addition, taking into account the property that is diagonalizable, the Markov transition probability matrix can be derived in closed form by the Lagrange-Sylvester interpolation (Sylvester, 1883) as where is the Frobenius covariant, which is independent of and defined as In the actual computation, whichever is more convenient can be used, either Equation (13) or Equation (15). Two important features of the model proposed in this study are worth highlighting here. First, the proposed model can be considered an improved form of the stagedhomogeneous Markov model originally used in Butt et al. (1987). Instead of a deterministically specified transition rate or transition probability, the proposed model represents the dynamic transition of the stage probabilistically using a latent CTMC. Second, the proposed model can be interpreted as a type of semi-Markov models because it satisfies the property of the semi-Markov process in which the entire stochastic process is not Markovian, but only at the instant of a jump (transition of the DC), it is Markovian. However, the commonly used semi-Markov models (Kobayashi et al., 2010;Thomas & Sobanjo, 2013) do not have an analytic solution for transition probabilities because the probability distribution of the holding time (sojourn time) is represented by a Weibull distribution. On the other hand, the proposed model has transition probabilities that are expressed analytically, and this is a major difference between the proposed model and existing semi-Markov models. Note that although the stochastic process becomes a CTMC when the probability distribution of the holding time in the semi-Markov model is exponential, this means that it becomes a CTMC without the LS, which is different from the proposed model.

Posterior probability density
Here, Bayesian estimation is used to obtain an estimate of parameter . When the joint prior probability density of is defined as ( ), the joint posterior probability density of F I G U R E 3 Parameter estimation flow of the proposed model with given observation can be given as where ( ) = ∫ ( | ) ( )d is the normalizing constant, and is the parameter space. Since it is difficult to obtain ( ) and thus ( | ) analytically, samples following the posterior probability density are numerically obtained with the MCMC method. Specifically, Gibbs sampling (Geman & Geman, 1984) is used, which samples each element of , that is, , 1 , 2 ( ∈ ; 1 , 2 ∈ , 1 < 2 ) or 1 , 2 ( 1 , 2 ∈ , 1 < 2 ), from the conditional posterior probability density. Since the conditional posterior probability density cannot be given analytically even if Gibbs sampling is used, the parameters are sampled from each conditional posterior probability density using the random-walk Metropolis-Hastings algorithm (Hastings, 1970). The number of samples of parameters to be obtained by the MCMC method is set by trial and error to a sufficient number such that it can be judged on the basis of Geweke test statistics (Geweke, 1996) that the obtained samples have converged to the posterior distribution. Note that the truncated normal distribution with support (0, ∞] is used as the proposal distribution of the random-walk Metropolis-Hastings algorithm because the space of each parameter is not ℝ but ℝ + . The parameter estimation flow of the proposed model is shown in Figure 3.

Condition setting
A numerical example is performed by re-estimating the parameter from artificially generated data with a given , to verify the feasibility of parameter estimation, in par-ticular the identifiability of the parameters in the proposed model. It is supposed that = 4, = 3, and the transition rate matrix is given as shown in Table 1. Ten DCs at 3-year intervals, which are randomly generated with the above , are considered as one sample, and the generated 5,000,000 samples are used as the database to estimate the parameters. Assuming noninformative priors, parameters 1,1,2 , 1,2,3 , 1,3,4 , 2,1,2 , 2,2,3 , 2,3,4 , 3,1,2 , 3,2,3 , 3,3,4 , 1,2 , and 2,3 are estimated from the database by Gibbs sampling from each conditional posterior using the random-walk Metropolis-Hastings algorithm. Five thousand samples for each parameter are obtained by the MCMC method, with the first 1000 samples removed as burn-in, and the remaining 4000 samples are considered as samples from the posterior probability density. The means of the sampled parameters are used as the parameter estimates.

Estimated result
The parameter estimates, 90% credible intervals, and Geweke test statistics are listed in Table 2. Since the absolute values of the Geweke test statistics for all parameters are less than 1.96, the hypothesis that the sampled parameters are convergent at the 5% significance level cannot be rejected. This leads us to judge that all parameters converge to the posterior distributions. For all 11 parameters, the estimates and the parameter values set in Table 1 are similar and all the true parameter values are within the 90% credible intervals. This implies that the parameters have been estimated correctly and the identifiability of the parameters in the proposed model is satisfied in the ideal situations where a sufficient sample is obtained. However, it should be noted that when the sample size is small, the identifiability may not be satisfied due to the small amount of data. In such cases, it may be necessary to exclude some of the parameters from the estimation.
In addition, models with different were estimated from the same data set, respectively, to examine whether the model with the correct number of LSs is estimated as the statistically most desirable model. The WAIC (Watanabe-Akaike information criterion) (Watanabe, 2010) for each model with different is shown in Figure 4. As can be seen from the figure, the WAIC is minimum for = 3, that is, the model with the correct set for data generation is estimated as the statistically most desirable model.

6
REAL-WORLD EXAMPLE

Condition setting
Here, an open data set (Ontario's Open Data Team, 2018) from bridge conditions observed in Ontario, Canada, is  used for a real-world example to validate the estimation capability of the proposed method for the time-dependent deterioration process of real infrastructure assets. In the data set, DCs of bridges have been recorded as BCI (bridge condition index), which takes a continuous value from 100 (good condition) to 0 (poor condition). The DC of a bridge at a time of inspection is represented by a single BCI value. The original data set includes historical BCI records, time of inspection, and time of construction for each bridge.
Here, it is assumed that BCI is 100 at the time of construction of the bridge. Statistics of the data, definitions of discrete state spaces, and so forth are shown in Table 3. Bridges that have been repaired or rehabilitated in the past are excluded from the data set.
Due to the small sample size, it is assumed that 1,1,2 = 1,1,3 = 2,1,3 = 0 and 2,2,3 = 1,2,3 and estimate only 2,1,2 , 1,2,3 , and 1,2 . Note that estimation of parameters assuming 3 or greater was also tried; however, the convergence of the parameters was not satisfactory. Assuming noninformative priors, parameters 2,1,2 , 1,2,3 , and 1,2 are estimated from the generated data by Gibbs sampling from each conditional posterior using the random-walk Metropolis-Hastings algorithm. Ten thousand samples for each parameter are obtained by the MCMC method, with the first 1000 samples removed as burn-in, and the remaining 9000 samples are considered as samples from the posterior distribution. The parameter means in the posterior distribution are used as the parameter estimates.
Here, the authors discuss the significance of modeling the deterioration process using the deterioration state (BCI) for the entire bridge in this real-world example, instead of modeling the deterioration process for each bridge element. Recent bridge management systems may often have a deterioration prediction process for each bridge element (El Hajj et al., 2017;Huang et al., 2010). This process is useful in developing detailed budget and repair plans for individual bridges. On the other hand, asset owners who manage a large number of bridges often have other processes within their asset management system to develop budget plans and evaluate performance, from a schematic, long-term perspective for the entire asset portfolio. In such processes, it cannot be denied that TA B L E 3 Overview of data for modeling of actual bridge deterioration predicting deterioration for each bridge element may result in an overly detailed analysis. In addition, it is suggested that prediction of entire bridge deterioration is still meaningful because (1) BCI is still being recorded in the Ontario data set used in this study and (ii) analyses that predicted deterioration by BCI for entire bridges have been actively conducted in recent years (Taghaddos & Mohamed, 2019). For these reasons, the authors decided to use the BCI for the entire bridge as a deterioration index in this real-world example, while also recognizing the importance of predicting deterioration for each bridge element. It should be noted that the proposed method can be used for bridge elements and other types of infrastructure assets as long as DC time series are available in the same format as in this real-world example.

Estimated result
The parameter estimates, 90% credible intervals, and Geweke test statistics are listed in Table 4. The mean value of the log-likelihood over 9000 samplings, excluding burnin, is also shown in the table. For the same reasons as in the numerical example, it can be judged that all parameters converge to the posterior distributions.
Using the parameter estimates, the Markov transition probabilities can be calculated as follows: Estimated Markov transition probability: Comparing Equation (24) and Equation (25), it can be confirmed that the transition probability of DC 1 takes a larger value, that is, the deterioration progresses faster, when LS is 2 than when LS is 1. In the estimation results, the Markov   Figure 5, states 2 and 3 are unreachable due to the assumptions set in the Markov transition probabilities. From these figures, it can be confirmed that a time-dependent deterioration process, in which the LS transitions from 1 to 2 in the middle of the service life of the bridge (within approximately 25 years of bridge age), that is, a deterioration stage shifted to a stage in which the deterioration of the DC progresses faster, has been inferred. Note that the shape of the deterioration process shown in Figure 6 can vary depending on the definition of the DC (in this example, the threshold value of BCI). For example, by changing the definition of the DC, more general infrastructure deterioration curves developed using human (inspector) ratings (reverse sigmoidal curve) might be estimated. Although the validity of the definition of the DC and its usefulness in asset management can be discussed through life cycle cost analysis or risk analysis using the estimated result of the deterioration process, it is out of the scope of this study.

Comparison
In this subsection, the traditional homogeneous CTMC and the semi-Markov model are taken as the existing and state-of-the-art methods, respectively, and compared with the proposed method. When  = {1} in the setting of the proposed method, the proposed method becomes consistent with the conventional time homogeneous CTMC. The result of estimating only the parameters 1,1,2 and 1,2,3 with  = {1} is used as the result of the existing method. The parameter estimates with the existing method are listed in Table 5, and the Markov transition probability of the DC is calculated with the parameter estimates with the existing method as The WAIC of the proposed method in Table 4 is lower than that of the existing method in Table 5, which indicates the statistical superiority of the proposed model. The time evolution of the DC state vector estimated with the existing method is shown in Figure 8. In the figure, the boundaries of each DC in Figure 6 are shown as dashed lines. From the figure, it can be seen that the transition processes of the estimated state vectors are obviously different between the proposed method and the existing method. In particular, the existing method underestimates the occurrence probabilities of DC 2 and DC 3 in the latter part of the bridge life (after approximately 20 years of bridge age), which can be attributed to the fact that the deterioration process estimated by the existing method is not time dependent. In Taghaddos and Mohamed (2019), a condition with a BCI value of 70 or less is defined as a condition requiring maintenance. In this study, although it is a slightly pessimistic assessment, it is supposed that DC 3, which has a BCI of 75 or less as shown in Table 3, is a condition that requires maintenance. This and the results of Figure 8 suggest that in this real-world example, from the bridge management perspective, the existing model risks overestimating the amount and cost of necessary repairs in the near future (from the present to approximately 10 years from now) while underestimating the amount and cost of necessary repairs in the more distant future (after 20 years).
The percentages of the DCs in the observed data for each year since construction are shown in Figure 9. In Figure 9a and b, the boundaries of each DC in Figure 6 and Figure 8 are shown as dashed lines, respectively. From these figures, it can be visually confirmed that the proposed model is able to estimate the deterioration process more similar to the observed data than the existing model (homogeneous CTMC).

Comparison with state of the art
The semi-Markov model is considered the state-of-the-art method, which is supposed to be widely used in current bridge management systems, and compared with the proposed method. The formulation and estimation method of the continuous-time semi-Markov model used in this study are given in Appendix B. The parameter estimates with the state-of-the-art method (i.e., semi-Markov model) using the same data set as 6.2 and 6.3.1 are listed in Table 6. Since the WAIC was lower when 2 = 1 was given exogenously as a deterministic value, a model was adopted in which 2 was not included as a parameter to be estimated. Although the mean value of the log-likelihood is slightly larger for the semi-Markov model than for the proposed model, the WAIC of the proposed method shown in Table 4 is lower than that of the state-of-the-art method shown in Table 6, which may suggest the statistical superiority of the proposed method. It should be noted, however, that in the estimation results of the semi-Markov model, the WAIC TA B L E 6 Parameter estimates of the real-world example: The state-of-the-art method The time evolution of the DC state vector estimated with the state-of-the-art method is shown in Figure 10. The time evolution of the DC state vector can be calculated analytically when using the proposed method, whereas simulation or numerical integration is required when using the state-of-the-art method. This figure shows Computation times required to obtain 10,000 samples from the posterior distribution for each parameter are shown in Table 7 for the proposed method and for the state-of-the-art method, respectively. All calculations in this paper were done in MATLAB R2019a on a computer with the following features: OS: Windows 10 Professional 64-bit; CPU: Intel Core i9 9900X 10core/20thread 3.5 GHz; Memory: 16.00 GB. MATLAB functions "integral" and "integral2" were used to numerically calculate the values of the likelihood function of the semi-Markov model. The proposed method reduces the computation time to (3, 633.605∕7, 500.421) × 100 = 48.45 (%) of the state-of-the-art method.

Forecasting future DC with the estimated result
This subsection explains how the estimation results of the proposed model can be used to forecast future deterioration processes based on observed information. Suppose that DCs have been observed at each of the times of inspection (calendar time) 1 , 2 , … , for an asset whose time of construction is 0 . Assume that 1 < 2 holds when 1 < 2 . The inspection interval is denoted by = − −1 . For timing ( = 0, … , ), denote the observed DC by and the unobserved LS by . Suppose that 0 = 1 and 0 = 1, and represents the current timing. The observation information set is defined as = { 0 , … , }.
The occurrence probability of LS ∈  at current timing with the given observations, ( | ), can be represented as the conditional occurrence probability of LS given , as ( , ), ( +1 , +1 ) ( +1 ) +1 ( +1 ) ( < − 1) where ( , ) is the joint occurrence probability of and , and ( ) is the occurrence probability of , which can be given by marginalizing ( , ) with respect to . Note that the Markov transition probabilities in the above equations are assumed to be calculated using the parameter estimates. The conditional occurrence probabilities of future state ∈  and future DC ∈  at + ( > 0) with given , ( | ) and ( | ), can be analytically derived, respectively, as Here, ( | ) is calculated based on the estimated Markov transition probabilities in a situation where = {1, 1}. Figure 11 shows the variation of the probability 10 (3 | {1, 1}) that the DC will be 3 in 10 years ( = 10), in response to the current bridge age . Note that the proposed model can also quantify differences in future deterioration processes due to the path dependence of past DCs (i.e., differences in ), in addition to bridge age. However, the path dependence is not considered here because the proposed model estimated in this real-world example has a model structure in which the LS becomes 2 deterministically if the DC transitions to 2. From Figure 11, it can be seen that differences in the probability of deterioration at a future time point due to differences in the current bridge age can be quantified by the proposed method, which is capable of representing time-dependent deterioration processes. On the other hand, the existing method estimates the probability of deterioration at a future time point as a constant value regardless of the bridge age. Figure 12 shows the variation of the probability (3 | {1, 1}) with respect to target future year for forecasting . For example, comparing "Proposed: Time of construction" and "Existing," it can be seen that the prediction results of the proposed method have a relatively smaller occurrence probability of DC 3 in the near future and a relatively larger occurrence probability of DC 3 in the far future than the prediction results of the existing method, which represents the deterioration process with a constant DC transition rate over time. This means that the proposed method is able to express the accelerated increase in future deterioration probability. Based on the above, it can be considered that the proposed method can forecast the future flexible deterioration F I G U R E 1 2 Variation of (3 | {1, 1}) in response to target future year process, while the existing method is likely to under-or overestimate the probability of future deterioration.

DISCUSSION
In this research, a time-inhomogeneous CTMC with analytic transition probabilities has been proposed. The analytic transition probabilities in the proposed model can (1) improve the mathematical tractability of the model and (2) reduce the computation time for statistical model estimation and future deterioration prediction. In addition, it has been confirmed that the identifiability of the model parameters is satisfied through the simulation example in Section 5. The effectiveness of the proposed model for existing studies has been demonstrated by the real-world example in Section 6. As shown in Section 6.3.1, it has been confirmed based on the comparison of the WAIC that the proposed time-inhomogeneous CTMC can estimate statistically superior deterioration prediction results in comparison with the existing homogeneous CTMC. As shown in Section 6.3.2, through comparison with the stateof-the-art inhomogeneous CTMC (semi-Markov model), it has been empirically confirmed that the proposed model can estimate a model with the equivalent accuracy to the state-of-the-art model in approximately 48% of the computation time of the state-of-the-art model.
On the other hand, there are two main limitations with respect to the proposed model. First, the model is constrained to yield the analytic transition probabilities. Future innovations in reducing numerical integration errors and improving computational efficiency could reduce the WAIC and computation time of the state-ofthe-art semi-Markov model in Section 6.3.2. Moreover, the possibility cannot be denied that there exists another undiscovered inhomogeneous CTMC with analytic transition probabilities whose statistical fit is better than that of the proposed model. Therefore, it is necessary to continuously study methods of modeling time-dependent deterioration processes, taking into consideration both computational simplicity and modeling accuracy. Second, the validity of the proposed model has been demonstrated only in one example. Since some of the findings obtained from the real-world example are applicable only to the bridge group analyzed, it is important to apply the proposed model to other bridge databases and to various types of assets to investigate its applicability.

CONCLUSION
In this research, it has been revealed that a regimeswitching model structure incorporating an unobservable LS enables the formulation of a time-inhomogeneous CTMC with analytic transition probabilities for representing time-dependent deterioration processes. The proposed model has been verified and validated through simulation and real-world examples, respectively. The proposed method can be used to estimate long-term asset performance and budgets with simplicity and accuracy. In addition, the analytic transition probabilities of the proposed model could be exploited to reduce the computational burden of optimizing intervention strategies. Therefore, a major possible future research is to develop a decisionmaking model based on the time-dependent deterioration process estimated by the proposed model, in which a partially observable Markov decision process based on ( | ) or ( | ) may be a good option.

A C K N O W L E D G M E N T S
This work was supported by the Japan Society for the Promotion of Science (JSPS) KAKENHI (grant numbers: JP19H00777 and JP20H02264). The data analyzed in Section 6 was provided by Ontario's Open Data Team (downloaded from https://data.ontario.ca/dataset/bridgeconditions), Government of Ontario, Canada, under the terms of the Open Government Licence -Ontario (OGL). : Transition rate matrix of the DC when the LS is , 1 , 2 : 1 × 2 element of ( 1 , 2 ∈ ) : Transition rate matrix of the LS 1 , 2 : 1 × 2 element of ( 1 , 2 ∈ )

R E F E R E N C E S
: State of the asset, which is defined as a state augmentation of the DC and LS : State space of the state : Transition rate matrix of the state : Diagonal element matrix of , , : Diagonal element of 1 , 2 : Diagonal matrix whose diagonal elements are all 1 , 2 : × identity matrix • The symbols defined in Section 3.2: ( ): Markov transition probability matrix of the state evolution 1 , 2 ( ): Markov transition probability from state 1 to 2 , which is 1 × 2 element of ( ) ( 1 , 2 ∈ ) : Time interval : Positive real number 1 , 2 : 1 × 2 element of : Diagonal matrix whose × element is , , : Eigenvector of , which corresponds to the eigenvalue , : × invertible matrix whose th column is , e, : Diagonal matrix whose × element is exp( , ) , 1 , 2 : 1 × 2 element of , 1 , 2 : 1 × 2 element of −1 : Frobenius covariant  Kobayashi et al. (2010) is called the "multistage Weibull hazard model" and not the "semi-Markov model" in the literature but can be interpreted as a continuous-time semi-Markov model. The details of the model are given in the literature, but, in the following, the model is formulated for the case where = 3.
Assuming transitions only to adjacent DCs, the probability density of the sojourn time in DC ∈  = {1, 2, 3} is expressed by the Weibull distribution as where > 0 is the shape parameter and > 0 is the scale parameter. The survival function is expressed as The joint occurrence probability of the elements of ℎ is defined as The likelihood can be given as Assuming noninformative priors, 10,000 samples for each parameter are obtained by the same MCMC method as described in Section 6.1. The integrals contained in Equation (B6) do not have analytic solutions and require approximate computation by numerical integration.