How to generate samples of ARIMA(p,d,q) model within an interval? I am want to generate samples from an ARIMA(p,d,q) or ARMA(p,q) model. There is a Python Package to generate ARMA samples. The problem is that I want to generate scenarios for demand which should be non-negative and also, usually there is an upper bound for demand. I would be thankful if can help me with how to generate samples from ARIMA(p,d,q) which fall in the interval [l,u].
 A: The first question is whether a model from the ARIMA class is appropriate for your data generating process. But I will leave it aside in this answer and only address your original question.

An ARIMA model will generally produce an anbounded sequence. One way to bound it would be to only use an MA model combined with a bounded error distribution (e.g. uniform, but not only; see "Distributions bounded on both sides" for more ideas). If your error term $\varepsilon$ is bounded as $\varepsilon\in[a,b]$ and your MA(q) model is
$$
y_t=\theta_0+\varepsilon_t+\theta_1\varepsilon_{t-1}+\dots+\theta_q\varepsilon_{t-q},
$$
then $y_t$ is bounded as follows:
$$
\theta_0-\xi 
\leq y_t\leq
\theta_0+\xi
$$
where $$
\xi=a\sum_{i=1}^q\min\{\theta_i,0\}+b\sum_{i=1}^q\max\{\theta_i,0\}.
$$
If you want $y_t$ to be bounded between $l$ and $u$, you need to solve
$\theta_0-\xi=l$ and $\theta_0+\xi=u$ for $a$ and $b$.

Another way of bounding ARIMA processes with AR terms and/or with unbounded error distributions would be to transform them, but that would no longer be an ARIMA process. The transformation could be done by using a cumulative density function (CDF) of a random variable with unbounded support. This would produce a number between $0$ and $1$. (Think logistic regression or probit model.) Then you could multiply the result by a scale parameter $u-l$ and add a location parameter $l$ to make them lie within $(u,l)$.
How would fit such a model in practice? You would take your raw data and do the opposite transformation. First, subtract $l$ and divide by $u-l$, then apply an inverse CDF to make them lie in $(\infty,+\infty)$ where they can be modelled with an arbitrary ARIMA. Then use the fitted model to do predictions or simulations and then transform these back to $[l,u]$ as described in the preceding paragraph.
