My line of thought goes like this. As per the AG model, each individual is represented by a counting process $N_i(t)$ with intensity $\lambda_i(t)$ which can be written as $\lambda_i(t) = \lambda_0(t) \exp(\beta' x_i(t))$. By this it is implicit that only the "current" value of $x_i$ matters for the intensity. The counting process grows when an individual has an event. I also assume independent censoring (actually stopping of the process). Note that these assumptions are implicit from the way that you fitted the model.
The probability of no events in $(a,b)$ (improperly, a "survival") is then
$$
S_i(t|s) = \exp \left( -\int_s^t \lambda_0(u) \exp(\beta'x_i(u)) du \right)
$$
Note that the probability of no events does not directly relate to the probability of one event, as the complementary event to "no events in a period" is "at least one event during a period".
Then say you are interested in $S(t|s_0)$. Then the idea would be to fix $x_i = x_i(s_0)$ (i.e. assume the covariates don't change after $s_0$), and we would get
$$
S_i(t|s_0) = S_i(t) / S_i(s_0)
$$
Since the part before $s_0$ cancels out, together with my assumption that only the current value of $x_i$ matters for the intensity, means that this is equal to
$$
S_i(t|s_i) = \exp \left( -\exp(\beta'x_i(s_0)) \int_{s_0}^t \lambda_0(u)du \right)
$$
Then in R you can do it like that (I give an example with a data set):
library(frailtypack)
data(readmission)
mod1 <- coxph(Surv(t.start, t.stop, event) ~ sex + charlson + cluster(id), data = readmission)
# set the covariate values at s_0
mycov <- data.frame(sex = "Female", charlson = "1-2")
sf <- survfit(mod1, newdata = mycov)
# with different s_0
par(mfrow=c(2,2))
time <- c(300, 500, 700, 1000)
for(i in time) {
pos <- which.max(sf$time[sf$time <= i])
S_s0 <- sf$surv[pos]
with(sf, plot(time[pos:length(time)], surv[pos:length(surv)] / S_s0), type = "l")
}
Here you get the plots of the survival curves, which correspond to the values
$(t, S(t|s_0))$ for $t\geq s_0$.
The other comments that I would give on this are the following: it is difficult to talk about the distribution of the next event. This is because in the AG formulation the time since previous event does not play any role. In other words, if you would like to take that into account, more complicated stochastic models should be used, where for example you include the "previous number of events" as a time-dependent covariate. This does complicate things a lot and the interpretation of the estimated quantities is most likely very difficult.
The second comment I have is about the nature of the time dependent covariates. Mostly, the AG works nicely with "external" covariates, such as air pollution, or something which is not directly measured on the subject ("external" of the recurrent event process). This is mostly because the first expression I wrote here is the probability of no events during $(a,b)$ relies on the assumption that the number of events in any given interval is Poisson distributed. This is true if the covariates are external. A discussion on this can be found in several textbooks, for example in Cook & Lawless at section 2.5. If your time-dependent covariates do depend on the recurrent event process, then it should be modeled jointly with the recurrent events process.
Parameterizations can be tricky, and get even more complicated when you move from the the logarithmic time scale, often used to illustrate a parametric model, to an absolute time scale. I find myself re-learning each time I try to answer a question like this.
I find these course notes to be a helpful, succinct summary. They cover several models of the general form
$$\log T = \alpha + \sigma W ,$$
where $\alpha$ is a "location" parameter, $\sigma$ is a "scale" parameter, and $W$ has a specified probability distribution. For a log-logistic model, $W$ is standard logistic, with random values given by rlogis() in R.
The trick is to find how to express the parameters of a log-logistic survival model in terms of $\alpha$ and $\sigma$ above. Those course notes show that the survival function $S(t)$ for a log-logistic model takes the form
$$ S(t) = \frac{1}{1+(\lambda t)^p},$$
where $\alpha = - \log \lambda$ and $\sigma = 1/p$. Those two relationships between parameters in the $\log T$ scale and the $T$ scale hold for many model types discussed in the course notes.
To check this and see what values survreg() returns for coefficients, generate some random log-logistic survival times, following your general approach. I did this for $\lambda = 2$ and $p = 3$, plugging into the formula for $\log T$, then exponentiating:
set.seed(123)
lldeathtimes <- exp(-log(2)+ (1/3)*rlogis(1000))
loglog.fit <- survreg(Surv(lldeathtimes) ~ 1, dist = "loglogistic")
loglog.fit$icoef
#(Intercept) Log(scale)
# -0.6985667 -1.1065919
You can then see how the reported (Intercept) and Log(scale) values correspond to the parameters of the model. Note that the (Intercept) is close to $-\log(2)$, or $\alpha$ in the formula for $\log T$ when $\lambda = 2$. For Log(scale), note that the corresponding scale value is close to $1/3$, or $\sigma$ in the formula for $\log T$ when $p=3$.
If you plot the synthetic survival data and overlay the intended log-logistic survival function in the parameterization I used:
plot(survfit(Surv(lldeathtimes)~1),xlim=c(0,4))
curve(1/(1+(x*2)^3),from=0,to=4,col="blue",add=TRUE)
you will see that they agree well.
So for a model of the above general form for $\log T$, (Intercept) is the "location" parameter $\alpha$ and scale is the "scale" parameter.
What gets confusing, particularly for Weibull models, is that the rweibull() parameterization differs from the used by survreg() with even a switching of names.* From the manual page for survreg.distributions:
survreg scale parameter maps to 1/shape, linear predictor to log(scale)
So you can't count on the word "scale" to mean the same thing, even in a single software package. Combine that with the different ways to parameterize the underlying distributions (what I used from the course notes for log-logistic differs from what you got from Wikipedia) and you have a recipe for massive confusion. I'd recommend avoiding the use of words "shape" or "scale" or "location" to describe the parameters and just use symbols in equations, to be most specific about just what you mean.
*I think that you have $\alpha$ and $\beta$ interchanged in your formulas for your Weibull model. To check, plot the synthetic data and the corresponding intended survival curve, in the way I used for the log-logistic model. Your exponential model seems OK.
Best Answer
According to the mean you give, you use the following parametrisation for the Weibull distribution: $$ \textrm{if }X\sim \textrm{Weibull}(\lambda, \alpha) \textrm{ then } f_X(x) = \lambda \alpha x^{\alpha - 1} \exp(-\lambda x^\alpha), $$ with $\lambda > 0$ a scale parameter, and $\alpha > 0$ a shape parameter.
dweibull() from R, as well as wikipedia, use another parametrisation. The conversion is as follows: $$ \textrm{shape} = \alpha \quad \textrm{and} \quad \textrm{scale} = \left(\frac{1}{\lambda} \right)^{\tfrac{1}{\alpha}}, $$ where $\textrm{shape}$ and $\textrm{scale}$ are those given in dweibull() and wikipeida.
Let $\mathbf{x}'\mathbf{\beta} = x_1\beta_1 + x_2\beta_2 + \dotsb$ be the linear predictor.
Assuming a proportional hazards structure and a $\textrm{Weibull}(\lambda, \alpha)$ distribution at baseline, the hazard rate is written \begin{align*} h(t) & = h_0(t) \exp(\mathbf{x}'\mathbf{\beta}) \\ & = \lambda \alpha t^{\alpha - 1} \exp(\mathbf{x}'\mathbf{\beta}). \end{align*} The corresponding pdf is $$ f(t) = \lambda \alpha t^{\alpha - 1} \exp(\mathbf{x}'\mathbf{\beta}) \exp \left( - \lambda t^\alpha \exp(\mathbf{x}'\mathbf{\beta}) \right). $$ That is, $T$ has a Weibull distribution with the same shape $\alpha$ but the scale parameter is changed from $\lambda$ to $\lambda \exp(\mathbf{x}'\mathbf{\beta})$: $$ T \sim \textrm{Weibull}(\lambda \exp(\mathbf{x}'\mathbf{\beta}), \alpha) $$ and we have $$ E[T] = \frac{\Gamma(1 + \tfrac{1}{\alpha})}{\left(\lambda\exp(\mathbf{x}'\mathbf{\beta})\right)^{\tfrac{1}{\alpha}}}. $$
An example without covariate: