Here we have a arrival process. The inter-arrival time follows a shifted negative exponential distribution. The density function of the distribution is:
$$f(t)=\lambda e^{-\lambda(t-\theta)},\quad\text{ where }t\ge\theta$$
How to derive the variance of the number of arrivals in time period of $T$?
Assuming that the inter-arrivals say $X_n$ ($n \geqslant 1$) are independent, you have a renewal process, see e.g. this course, or the classical references quoted in it: the book by D.R. Cox Renewal Theory or the one by S. Karlin and H.M. Taylor A First Course in Stochastic Processes, vol. 1 chap. 5.
The $n$-th arrival time from $t=0$ is the sum $S_n:= X_1 + X_2 + \dots + X_n$ for a specific initial condition: when $t=0$ is an arrival time. Then $X_1$ is distributed as are the $X_n$ for $n > 1$. A variant takes a specific stationary distribution for the first arrival $X_1$, leading to the stationary renewal process. The initial condition yet has no impact in the long run.
Let $N(t)$ be the number of arrivals on $(0,\,t)$. When $t$ is large, a renewal theorem states that $N(t)$ is approximately normal with mean $t/\mu$ and variance $t \sigma^2 / \mu^3$ where $\mu$ and $\sigma^2$ are the inter-arrival mean and variance. In your case, the theorem applies with $\mu = \theta + 1/\lambda$ and $\sigma = 1/\lambda$.
The distribution of $N(t)$ can also be found by noticing that $\text{Pr}\{N(t) \geq n\} = F_n(t)$ where $F_n$ is the distribution function of the sum $S_n$, and thus $\text{Pr}\{N(t) = n\} = F_n(t) - F_{n+1}(t)$. In your case, $X_n = X_n^\star +\theta$ where $X_n^\star$ follows a standard exponential with mean $1/\lambda$, so $F_n(t) = F_n^\star(t-n\theta)$ where $F_n^\star$ is the distribution function of the sum $S_n^\star$ relative to the $X_k^\star$, and an explicit formula based on Erlang's distribution can be used in numerical computations. Assume an arrival at $t=0$ and let $t^\star := t -n \theta$; if $t^\star > 0$, then $$ F_n^\star(t^\star) = 1 - \sum_{k=0}^{n-1} e^{-\lambda t^\star} \frac{(\lambda t^\star)^k}{k!} = \text{Pr}\left\{N^\star \ge n \right\} $$ where $N^\star$ is Poisson with mean $\lambda t^\star$. A similar formula can be used for $F_{n+1}(t)$. The number of probability masses $\text{Pr}\{N(t) = n\}$ to be computed must be such that the total mass is close to $1$.
theta <- 0.4; lambda <- 1.0;
mu <- theta + 1 / lambda; sigma <- 1 / lambda
t <- 10;
## asymptotic 'Exp'ectation and 'Var'iance from the central limit renewal thm
aExpN <- t / mu
aVarN <- t * sigma^2 / mu^3
## compute the distribution: 'nMax' should be chosen suitably.
## Pr{ N(t) = n } is 'prob[n + 1]' since array indices are >= 1
nMax <- 100; prob <- rep(0, nMax + 1)
for (n in 0:nMax){
tStar <- t - n * theta
if (tStar > 0) {
prob[n + 1] <- prob[n + 1] +
ppois(n - 1, lambda = lambda * tStar, lower.tail = FALSE)
}
tStar <- t - (n + 1) * theta
if (tStar > 0) {
prob[n + 1] <- prob[n + 1] -
ppois(n, lambda = lambda * tStar, lower.tail = FALSE)
}
}
names(prob) <- 0:nMax
ExpN <- sum((0:nMax) * prob)
VarN <- sum((0:nMax)^2 * prob) - ExpN^2
## compute (estimate) expectation and variance using a simulation
nSim <- 500000
set.seed(12345) ## to be reproducible
X <- theta + matrix(rexp(100 * nSim, rate = lambda), nrow = nSim, ncol = 100)
Nsim <- apply(X, MARGIN = 1, FUN = function(x) { sum(cumsum(x) < t) } )
## compare empirical and numerical distributions
prob1 <- table(Nsim) / length(Nsim)
prob2 <- cbind(prob1, prob[names(prob1)])
colnames(prob2) <- c("sim", "num")
barplot(t(prob2), beside = TRUE, legend = TRUE,
main = sprintf(paste("distr. of the number of arrivals",
"lambda = %5.2f, theta = %5.2f"), lambda, theta))
## compare Expectation and variance
res <- rbind(asympt = c(aExpN, aVarN),
sim = c(mean(Nsim), var(Nsim)),
num = c(ExpN, VarN))
colnames(res) <- c("Exp", "Var")
res
