A discreet-time Markov process for which a transition probability matrix $P$ is independent of time can be represented, or approximated, with a continuous-time Markov chain (CMTC) with constant transition rates using exponential probability distributions. How do we obtain a continuous-time representation when the transition rates are not constant but depend on time a given by Weibull distribution?
For constant rates, we have for CMTC
$$P(t) = \exp(Qt)$$
(as a result of solving the forward equation) where $Q$ is the transition rate matrix with entries related to the exponential distribution parameters.
The above expression for $P(t)$ would probably not longer be valid if the rates are not constant in time. Is there a known solution to this "generalized" problem (since exponential is a special case of Weibull for k=1)?
Add 1
The forward equation for $P(t)$ is given by
$$ \frac{\partial P(t)}{\partial t} = P(t)Q(t). $$
If I understand it correctly, then the transition rate matrix $Q$ contains on the off-diagonal hazard rates of the competing underlying processes representing CTMC. If these processes follow exponential distribution, then these rates are constant and $Q$ is independent of time. However, if these hazard rates follow Weibull distribution then, I suppose, the off-diagonal elements of $Q$ are given by
$$ q_{i,j} = \beta_{i,j} k_{i,j} t^{k_{i,j}-1}. $$
It seems to me that with these form of the hazard rates, it is still posible to solve analytically the forward equation to get
$$ P(t) = \exp(A(t)) $$
where
$$ (A)_{i,j} = \beta_{i,j} t^{k_{i,j}}. $$
I don't really know if this is correct or how to go from here to find, for example, limiting distribution.
Add 2
The system I am considering has 4 states: 1 "alive" and 3 "dead", so the transition probability matrix looks something like
$$ \begin{pmatrix} p_{0,0}(t) & p_{0,1}(t) & p_{0,2}(t) & p_{0,3}(t) \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ \end{pmatrix} $$
so, if I am not mistaken, we have the following system of ODE's (assuming that the forward equation still holds in this case):
$$ \begin{bmatrix} p^\prime_{0,0} = p_{0,0} \cdot q_{0,0}(t) \\ p^\prime_{0,1} = p_{0,1} \cdot q_{0,1}(t) \\ p^\prime_{0,2} = p_{0,2} \cdot q_{0,2}(t) \\ p^\prime_{0,3} = p_{0,3} \cdot q_{0,3}(t) \\ \end{bmatrix} $$
Is this correct? Would this system be solvable for $q_{0,i}(t) : i>0 $ given by Weibull hazard rates?
Add 3
If the above is correct, then with $q_{0,i}(t) = \beta_i k_i t^{k_i - 1} : i>0$, I get that
$$ \frac{d p_{0,0}(t)}{dt} = p_{0,0}(t) \left(- \sum_{i=1}^3 \beta_i k_i t^{k_i - 1} \right) $$
which, with $p_{0,0}(t=0) = 1$, gives
$$ p_{0,0}(t) = \exp \left(- \sum_{i=1}^3 \beta_i t^{k_i} \right) $$
and
$$ p_{0,i}(t) = \int_0^t \beta_i k_i s^{k_i - 1} \exp \left(- \sum_{i=1}^3 \beta_i s^{k_i} \right) ds + c_i $$
If we had $k_i = k \; \forall \; i$, then, with $p_{0,i}(t=0) = 0 : i>0$, the above can be integrated to get
$$ p_{0,i}(t) = \frac{\beta_i}{\sum_{i=1}^3 \beta_i} - \frac{\beta_i}{\sum_{i=1}^3 \beta_i} \exp \left(- t^k \sum_{i=1}^3 \beta_i \right) $$
but, alas, the $k_i$s are not the same.
Add 4
I mostly interested in the probabilities for $t \le 1$ and I have $\beta_i < 1 \; \forall \; i$, so I am thinking of approximating the exponential $\exp \left(- \sum_{i=1}^3 \beta_i t^{k_i} \right)$ with a second-degree polynomial in $\sum_{i=1}^3 \beta_i t^{k_i}$. I am thinking about using something like shifted Legendre polynomials on [0, 1], but my upper limit may be less than 1 (actually, more than $-\sum_i \beta_i$), so might need to look into adjusting the coefficients for this. I am also not sure if calculating the coefficients for the expansion of $\exp(x)$ and substituting $x=\sum_{i=1}^3 \beta_i t^{k_i}$ will give the same approximation as calculating the coefficients on $\exp \left(- \sum_{i=1}^3 \beta_i t^{k_i} \right)$ directly.
In the scenario with
$$Q=\begin{bmatrix} -(r_1(t)+r_2(t)+r_3(t)) & r_1(t) & r_2(t) & r_3(t) \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{bmatrix}$$
you can solve the forward equation mostly explicitly. Assuming you start with $p_0=1$ (it is easy enough to adjust if this isn't the case), you have $p_0(t)=\exp \left ( -\int_0^t \sum_{i=1}^3 r_i(s) \, ds \right )$. For the others you have $p_i(t)=\int_0^t r_i(s) p_0(s) \, ds=\int_0^t r_i(s) \exp \left ( -\int_0^s \sum_{j=1}^3 r_j(u) \, du \right ) \, ds$. In your situation this results in some integrals that will need numerical evaluation, even for getting the equilibrium probabilities, except under some extremely restrictive assumptions about the exponents $k_{0,j}$. But this is a pretty inexpensive class of integrals to evaluate, since the $r_j$ can be analytically integrated themselves.
That said, you can make a little bit of analytical progress in your situation by changing variables in the outer integral to $v:=\int_0^s r_i(u) \, du=\beta_{0,i} s^{k_{0,i}}$. Then you have $\int_0^s r_j(u) du = \beta_{0,j} s^{k_{0,j}}=\beta_{0,j} \beta_{0,i}^{-\frac{k_{0,j}}{k_{0,i}}} v^{ \frac{k_{0,j}}{k_{0,i}}}$.
Therefore you end up with
$$p_i(t)=\int_0^{t'} \exp \left ( -v -\sum_{j \neq i} \beta_{0,j} \beta_{0,i}^{-\frac{k_{0,j}}{k_{0,i}}} v^{ \frac{k_{0,j}}{k_{0,i}}} \right ) dv$$
where $t'=\int_0^t r_i(s) ds$.
You still can't generally evaluate that analytically, but it does standardize the situation in some sense, especially if you are seeking the stationary distribution (so that $t'$ is just $\infty$).