Can a Non-Stationary ARIMA Series Be Simulated Using arima.sim in R? I am performing an experiment using classes of ARIMA models which I have to simulate for instance AR, MA, ARMA and possibly ARIMA using arima.sim() function in R. I am able to AR, MA, ARMA but ARIMA seems impossible as the coefficient of AR must be greater than 1 ($|\phi| > 1$) for the series to be non-stationary.
When $|\varphi| < 1$
Simulate AR(1) of 10 sample size with a coefficient of 0.8
# simulate ARIMA(1, 0, 0)
n <- 10
phi <- 0.8
set.seed(837530)
wn <- rnorm(n, mean = 0, sd = 1)
ar1 <- sqrt((wn[1])^2/(1-phi^2))
for(i in 2:n){
  ar1[i] <- ar1[i - 1] * phi + wn[I]
}

Here the arithmetic of this expression $\varphi = \sqrt((\varepsilon _{1})^2/(1-\varphi^2))$ is valid in the sense that $\varphi = \sqrt(0.1309045/0.36)$ which is $\sqrt0.3636236$
When $|\varphi| > 1$
Simulate AR(1) of 10 sample size with a coefficient of 0.8
# simulate ARIMA(1, 0, 0)
n <- 10
phi <- 1.2
set.seed(837530)
wn <- rnorm(n, mean = 0, sd = 1)
ar1 <- sqrt((wn[1])^2/(1-phi^2))

# Warning message:
#In sqrt((wn[1])^2/(1 - phi^2)) : NaNs produced

for(i in 2:n){
  ar1[i] <- ar1[i - 1] * phi + wn[I]
}
ar1
 #[1] NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN

Here the arithmetic of this expression $\varphi = \sqrt((\varepsilon _{1})^2/(1-\varphi^2))$ is not valid in the sense that $\varphi = \sqrt(0.1309045/-0.44)$ which is $\sqrt-0.2975102$
Similar experiment is performed using arima.sim() function
When $|\varphi| < 1$
set.seed(837530)
ma1 <- arima.sim(n = 10, model = list(ar = c(0.8), order = c(1, 0, 0)), sd = 1)

No error message
When $|\varphi| > 1$
set.seed(837530)
ma1 <- arima.sim(n = 10, model = list(ar = c(1.2), order = c(1, 0, 0)), sd = 1)
#Error in arima.sim(n = 10, model = list(ar = c(1.2), order = c(1, 0, 0)),  : 
  #'ar' part of model is not stationary

What I want to know
Can it be said that a non-stationary series can not be simulated using arima.sim?
 A: You can't use stats::arima.sim to simulate from an ARIMA process where the AR part implies nonstationary behavior; it gives an error that explicitly states that this is not supported by this function. You can have unit roots (specified in the difference part, not the AR part) however, so there exist nonstationary ARIMA processes which it can simulate.
arima.sim does not prevent this situation because of what you mention concerning $\frac{\sigma^2}{1-\phi^2}$ leading to a negative variance: it does not start from a draw from the distribution $X_0 \sim \mathcal{N}\left(0,\frac{\sigma^2}{1-\phi^2}\right)$  by default, it starts from $X_0 = 0$ instead, so this wouldn't matter.
The only reason I see for why arima.sim would choose not to support such ARIMA models is that it uses a rule of thumb to automatically determine an appropriate burn-in period at the beginning of simulation:
n.start <- p + q + ifelse(p > 0, ceiling(6/log(minroots)),  0)

minroots is the smallest norm of the roots of the characteristic polynomial. In the AR(1) case, this is just $|1/\phi|$, which would lead to a negative burn-in if it's less than 1. It of course doesn't make sense to even compute a burn-in period in a situation where $|\phi| > 1$, since the purpose of such a thing is to take a process that isn't stationary but that converges to one asymptotically, and crudely discard the first few steps so that what's left is closer in behavior to the stationary distribution. arima.sim could probably simply set the burn-in to zero in that case, but it doesn't. I wouldn't expect much support to change this behavior considering that what you're trying to do is 1) trivial to do by hand, 2) unlikely to be widely useful to other people and 3) a change to a function that is part of a core package.
The "meaning" of drawing $X_0 \sim \mathcal{N}\left(0,\frac{\sigma^2}{1-\phi^2}\right)$ is that, when $|\phi|<1$, this is the stationary distribution of the process: if you start from this distribution, then $X_t$ will have this same distribution at all $t$. If $|\phi|<1$ and you start from another distribution, there is a long-run distribution which you reach asymptotically and which is identical to the stationary distribution (you could use a burn-in period to sample it approximately). However, when $|\phi| > 1$, there is no stationary distribution, the process does not converge to any such distribution in the long-run, and there is no particular reason to start from any distribution; you have to specify it as part of your model.
A: You can simulate stationary ARMA models using the rGARMA function in the ts.extend package.  If you want to extend this to ARIMA models then all you have to do is to simulate the ARMA model and then add the required number of differencing steps.  Extensions to non-stationary time-series processes with explosive roots can be done, but it requires specification of starting points for the series.
