Let $\mathbf \varepsilon = (\varepsilon_1, \varepsilon_2, \ldots, \varepsilon_n)$, $\mathbf y = (y_1,y_2,\ldots, y_n)$, and $\mathbf{\mu} = (\mu, \mu, \ldots, \mu)$. The relationship among these vectors given in the question is
$$\mathbf y - \mathbf\mu = \mathbf{X}\mathbf \varepsilon$$
The matrix $\mathbf{X} = \mathbf{I} + \theta\mathbf{J}$ where $\mathbf{I}$ is the $n$ by $n$ identify matrix and $\mathbf J$ is the $n$ by $n$ matrix with $1$ on the immediate subdiagonal (that is, where $i = j+1$) and $0$ elsewhere.
The problem is that when $|\theta| \gt 1$, $\mathbf{X}$ can be practically singular. (Since $\mathbf{X}^{-1} = \mathbf{I} + \sum_{i=1}^{n-1} (-1)^i \theta^i \mathbf{J}^i$ and the $i^\text{th}$ powers of $\mathbf J$ are also subdiagonal matrices with $1$s (on the diagonal $i$ steps below the main one), the largest entry in the inverse of $\mathbf X$ is $(-1)^{n-1}\theta^{n-1}$. This can easily overflow floating point representations. It will obliterate all precision once $(n-1)\log_{2}|\theta| \gt 52$, which for $\theta=5$ occurs once $n\ge 24$.)
R will even tell you about this problem:
n <- 30
theta <- 5; mu <- 10; sigma <- 3
set.seed(17); epsilon <- rnorm(n, sd=sigma)
x <- diag(1, n) + matrix(c(0, diag(1, n)[-n^2]), n, n) * theta
solve(x, rep(0, n))
Error in solve.default(x, rep(0, n)) :
system is computationally singular: reciprocal condition number = 7.15828e-22
All statistical software provides a way around this: compute the generalized inverse (aka "least squares fit"). In R this can be done with lm (among other tools):
y <- mu + x %*% epsilon # Compute the series `y`
d <- as.data.frame(cbind(y-mu, x))# Create a data frame for computing the solution
fit <- lm(V1 ~ . - 1, data=d) # Obtain the solution
epsilon.hat <- coef(fit) # Extract it from `fit`
This time there are no complaints. Just for fun, I ran this procedure for $n=3000$. Here is a plot of the recovered values of $\mathbf \varepsilon$ against the original ones, with the diagonal line (of perfect equality) superimposed for reference:
plot(epsilon, epsilon.hat, col="#00000060")
abline(c(0,1), col="#e0000080", lwd=2)
The recovery is not perfect, as evidenced by the single off-diagonal plotting symbol. Experimentation suggests this is unavoidable; it is rare for the recovered errors exactly to equal the original ones. The problems occur at the very end of the sequence, as evidenced by the vertical departures from zero (marked by the horizontal gray line):
res <- epsilon.hat - epsilon # Residuals
plot(n-0:9, head(rev(res), 10), type="b",
xlab="Index (t)", ylab="Residual", main="Last residuals")
abline(h=0, col="Gray", lty=2)
The very last one is NA (because lm recognized problems and dropped a variable, I believe).
Nevertheless, for such an unstable series that's a pretty good showing.
