Variance of arrival process with shifted exponential distribution 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$?
 A: 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


