Repository: statsmodels/statsmodels

Issue:
`ARIMAResults.simulate` doesn't respect `random_state`
#### Describe the bug

The method `ARIMAResults.simulate` doesn't seem to respect the parameter `random_state`. It returns different results for a fixed random state parameter.

#### Code Sample


```python
from statsmodels.tsa.arima.model import ARIMA
import numpy as np

x = np.random.randn(100)
mod = ARIMA(x, order=(1, 0, 0))
res = mod.fit()

for i in range(5):
    print(i, res.simulate(1, random_state=7))
```


#### Expected Output

For me, this prints a different number for each of the iterations. My expectation, however, is that the same number is generated on each invocation of `simulate` because `random_state` is fixed.


#### Output of ``import statsmodels.api as sm; sm.show_versions()``

<details>

INSTALLED VERSIONS
------------------
Python: 3.10.12.final.0
OS: Darwin 23.3.0 Darwin Kernel Version 23.3.0: Wed Dec 20 21:30:44 PST 2023; root:xnu-10002.81.5~7/RELEASE_ARM64_T6000 arm64
byteorder: little
LC_ALL: en_US.UTF-8
LANG: en_US.UTF-8

statsmodels
===========

Installed: 0.14.1 (/Users/adaitche/homebrew/lib/python3.10/site-packages/statsmodels)

Required Dependencies
=====================

cython: Not installed
numpy: 1.24.1 (/Users/adaitche/homebrew/lib/python3.10/site-packages/numpy)
scipy: 1.10.0 (/Users/adaitche/homebrew/lib/python3.10/site-packages/scipy)
pandas: 1.5.3 (/Users/adaitche/homebrew/lib/python3.10/site-packages/pandas)
    dateutil: 2.8.2 (/Users/adaitche/homebrew/lib/python3.10/site-packages/dateutil)
patsy: 0.5.6 (/Users/adaitche/homebrew/lib/python3.10/site-packages/patsy)

Optional Dependencies
=====================

matplotlib: 3.6.3 (/Users/adaitche/homebrew/lib/python3.10/site-packages/matplotlib)
    backend: MacOSX
cvxopt: Not installed
joblib: Not installed

Developer Tools
================

IPython: 8.9.0 (/Users/adaitche/homebrew/lib/python3.10/site-packages/IPython)
    jinja2: 3.1.2 (/Users/adaitche/homebrew/lib/python3.10/site-packages/jinja2)
sphinx: Not installed
    pygments: 2.14.0 (/Users/adaitche/homebrew/lib/python3.10/site-packages/pygments)
pytest: 7.4.3 (/Users/adaitche/homebrew/lib/python3.10/site-packages/pytest)
virtualenv: Not installed

</details>


Discussion:
ChadFulton: Thanks for fixing this @bashtage!

Patch:
diff --git a/statsmodels/tsa/statespace/mlemodel.py b/statsmodels/tsa/statespace/mlemodel.py
index aa6f67a7e..08ad4c103 100644
--- a/statsmodels/tsa/statespace/mlemodel.py
+++ b/statsmodels/tsa/statespace/mlemodel.py
@@ -3644,6 +3644,10 @@ class MLEResults(tsbase.TimeSeriesModelResults):
         if iloc > self.nobs:
             raise ValueError('Cannot anchor simulation outside of the sample.')
 
+        # GH 9162
+        from statsmodels.tsa.statespace import simulation_smoother
+        random_state = simulation_smoother.check_random_state(random_state)
+
         # Setup the initial state
         if initial_state is None:
             initial_state_moments = (
@@ -3652,7 +3656,7 @@ class MLEResults(tsbase.TimeSeriesModelResults):
 
             _repetitions = 1 if repetitions is None else repetitions
 
-            initial_state = np.random.multivariate_normal(
+            initial_state = random_state.multivariate_normal(
                 *initial_state_moments, size=_repetitions).T
 
         scale = self.scale if self.filter_results.filter_concentrated else None


Test patch:
diff --git a/statsmodels/tsa/arima/tests/test_model.py b/statsmodels/tsa/arima/tests/test_model.py
index d04a5c36c..fb0502b26 100644
--- a/statsmodels/tsa/arima/tests/test_model.py
+++ b/statsmodels/tsa/arima/tests/test_model.py
@@ -413,3 +413,22 @@ def test_hannan_rissanen_with_fixed_params(ar_order, ma_order, fixed_params):
         res = mod.fit(method='hannan_rissanen')
 
     assert_allclose(res.params, desired_p.params)
+
+
+@pytest.mark.parametrize(
+    "random_state_type", [7, np.random.RandomState, np.random.default_rng]
+)
+def test_reproducible_simulation(random_state_type):
+    x = np.random.randn(100)
+    res = ARIMA(x, order=(1, 0, 0)).fit()
+
+    def get_random_state(val):
+        if isinstance(random_state_type, int):
+            return 7
+        return random_state_type(7)
+
+    random_state = get_random_state(random_state_type)
+    sim1 = res.simulate(1, random_state=random_state)
+    random_state = get_random_state(random_state_type)
+    sim2 = res.simulate(1, random_state=random_state)
+    assert_allclose(sim1, sim2)