Solution (source code)

= Solution

For transition-intensity vector $q$ and generator $Q(q)$, the expected infectious occupancy over three months from $S$ is the <Markov reward model> quantity
$$
m(q)=\mathbb E_S\int_0^3\mathbf1_{\{X_t=I\}}dt
=\int_0^3[e^{Q(q)t}]_{SI}\,dt.
$$
Insert the fitted intensities to obtain $m(\widehat q)$, evaluating the matrix exponential and integral numerically if necessary.

If $\eta=(\log q_{SE},\log q_{EI},\log q_{IR})$ has estimated covariance $V$, calculate the numerical gradient $g=\nabla_\eta m(e^\eta)$ at $\widehat\eta$. The <delta method> gives estimated variance $g^TVg$ and the approximate confidence interval
$$
m(\widehat q)\pm z_{1-\alpha/2}\sqrt{g^TVg}.
$$
Simulation from $N(\widehat\eta,V)$ followed by transformation through $m$ gives a useful asymmetric alternative.