diff --git a/development/modules/auxiliary_reliability.py b/development/modules/auxiliary_reliability.py
index 1db71d38..19d21543 100644
--- a/development/modules/auxiliary_reliability.py
+++ b/development/modules/auxiliary_reliability.py
@@ -3,7 +3,6 @@ import shlex
 from copy import deepcopy
 
 import numpy as np
-from statsmodels.tools.eval_measures import rmse
 
 import respy
 from development.modules.auxiliary_shared import aggregate_information
@@ -11,6 +10,7 @@ from development.modules.auxiliary_shared import cleanup
 from development.modules.auxiliary_shared import send_notification
 from development.modules.auxiliary_shared import update_class_instance
 from development.modules.config import SPEC_DIR
+from respy.python.shared.shared_auxiliary import rmse
 from respy.scripts.scripts_compare import scripts_compare
 from respy.scripts.scripts_simulate import scripts_simulate
 
diff --git a/environment.yml b/environment.yml
index e041d1c8..d1be3669 100644
--- a/environment.yml
+++ b/environment.yml
@@ -8,7 +8,6 @@ dependencies:
   - pandas=0.24
   - numba=0.43
   - scipy=1.2.1
-  - statsmodels=0.9.0
   - pip:
     - codecov
     - doc8
diff --git a/respy/python/shared/shared_auxiliary.py b/respy/python/shared/shared_auxiliary.py
index 1e6eeef6..6ae9e8a9 100644
--- a/respy/python/shared/shared_auxiliary.py
+++ b/respy/python/shared/shared_auxiliary.py
@@ -923,3 +923,78 @@ def number_of_triangular_elements_to_dimensio(num):
 
     """
     return int(np.sqrt(8 * num + 1) / 2 - 0.5)
+
+
+def ols(y, x):
+    """Ols implementation using a pseudo inverse.
+
+    Parameters
+    ----------
+    x (ndarray): n x n matrix of independent variables.
+    y (array): n x 1 matrix with dependant variable.
+
+    Returns
+    -------
+        beta (array): n x 1 array of estimated parameter vector
+
+    """
+    beta = np.dot(np.linalg.pinv(x.T.dot(x)), x.T.dot(y))
+    return beta
+
+
+def mse(x1, x2, axis=0):
+    """mean squared error.
+
+    Parameters
+    ----------
+    x1, x2 : array_like
+       The performance measure depends on the difference between these two
+       arrays.
+    axis : int
+       axis along which the summary statistic is calculated
+
+    Returns
+    -------
+    mse : ndarray or float
+       mean squared error along given axis.
+
+    Notes
+    -----
+    If ``x1`` and ``x2`` have different shapes, then they need to broadcast.
+    This uses ``numpy.asanyarray`` to convert the input. Whether this is the
+    desired result or not depends on the array subclass, for example
+    numpy matrices will silently produce an incorrect result.
+
+    """
+    x1 = np.asanyarray(x1)
+    x2 = np.asanyarray(x2)
+    return np.mean((x1 - x2) ** 2, axis=axis)
+
+
+def rmse(x1, x2, axis=0):
+    """root mean squared error
+
+    Parameters
+    ----------
+    x1, x2 : array_like
+       The performance measure depends on the difference between these two
+       arrays.
+    axis : int
+       axis along which the summary statistic is calculated
+
+    Returns
+    -------
+    rmse : ndarray or float
+       root mean squared error along given axis.
+
+    Notes
+    -----
+    If ``x1`` and ``x2`` have different shapes, then they need to broadcast.
+    This uses ``numpy.asanyarray`` to convert the input. Whether this is the
+    desired result or not depends on the array subclass, for example
+    numpy matrices will silently produce an incorrect result.
+
+    """
+    x1 = np.asanyarray(x1)
+    x2 = np.asanyarray(x2)
+    return np.sqrt(mse(x1, x2, axis=axis))
diff --git a/respy/python/solve/solve_auxiliary.py b/respy/python/solve/solve_auxiliary.py
index 9e8d45fd..4c30800a 100644
--- a/respy/python/solve/solve_auxiliary.py
+++ b/respy/python/solve/solve_auxiliary.py
@@ -3,17 +3,16 @@ import shlex
 
 import numpy as np
 import pandas as pd
-import statsmodels.api as sm
 from numba import njit
 
 from respy.custom_exceptions import StateSpaceError
-from respy.python.record.record_solution import record_solution_prediction
 from respy.python.record.record_solution import record_solution_progress
 from respy.python.shared.shared_auxiliary import calculate_rewards_common
 from respy.python.shared.shared_auxiliary import calculate_rewards_general
 from respy.python.shared.shared_auxiliary import create_covariates
 from respy.python.shared.shared_auxiliary import get_continuation_value
 from respy.python.shared.shared_auxiliary import get_emaxs_of_subsequent_period
+from respy.python.shared.shared_auxiliary import ols
 from respy.python.shared.shared_auxiliary import transform_disturbances
 from respy.python.shared.shared_constants import HUGE_FLOAT
 from respy.python.shared.shared_constants import MISSING_FLOAT
@@ -398,9 +397,7 @@ def pyth_backward_induction(
             # Create prediction model based on the random subset of points where the
             # EMAX is actually simulated and thus dependent and independent variables
             # are available. For the interpolation points, the actual values are used.
-            emax = get_predictions(
-                endogenous, exogenous, max_emax, is_simulated, file_sim, is_write
-            )
+            emax = get_predictions(endogenous, exogenous, max_emax, is_simulated)
 
         else:
             emax = construct_emax_risk(
@@ -529,6 +526,7 @@ def get_endogenous_variable(
         state has reached maximum education.
 
     """
+
     emax = construct_emax_risk(
         rewards[:, -2:], rewards[:, :4], emaxs, draws_emax_risk, delta, max_education
     )
@@ -538,53 +536,46 @@ def get_endogenous_variable(
     return endogenous
 
 
-def get_predictions(endogenous, exogenous, max_emax, is_simulated, file_sim, is_write):
-    """Fit an OLS regression of the exogenous variables on the endogenous variables and
-    use the results to predict the endogenous variables for all points in the state
-    space.
+def get_predictions(endogenous, exogenous, maxe, is_simulated):
+    """ Get ols predictions.
+
+    Fit an OLS regression of the exogenous variables on the endogenous variables and
+    use the results to predict the endogenous variables for all points in state space.
 
     Parameters
     ----------
-    endogenous : np.ndarray
-    exogenous : np.ndarray
-    max_emax : np.ndarray
-    is_simulated : np.ndarray
-    file_sim : ???
-    is_write : bool
+    endogenous (np.array): 1d numpy array with endogenous variable
+    exogenous (np.array): 2d numpy array with exogenous variables
+    maxe (np.array): 1d numpy array with the maximum expected value of the shocks,
+    has the same length as endogenous.
+    is_simulated (np.array): boolean array that is True at points for which we
+    calculate the exact solution.
 
     """
     # Define ordinary least squares model and fit to the data.
-    model = sm.OLS(endogenous[is_simulated], exogenous[is_simulated])
-    results = model.fit()
+    beta = ols(endogenous[is_simulated], exogenous[is_simulated])
 
-    # Use the model to predict EMAX for all states. As in Keane & Wolpin (1994),
-    # negative predictions are truncated to zero.
-    endogenous_predicted = results.predict(exogenous)
+    # Use the model to predict EMAX for all states. As in Keane & Wolpin (1994),negative
+    # predictions are truncated to zero.
+    endogenous_predicted = exogenous.dot(beta)
     endogenous_predicted = np.clip(endogenous_predicted, 0.00, None)
 
     # Construct predicted EMAX for all states and the replace interpolation points with
     # simulated values.
-    predictions = endogenous_predicted + max_emax
-
-    predictions[is_simulated] = endogenous[is_simulated] + max_emax[is_simulated]
+    predictions = endogenous_predicted + maxe
+    predictions[is_simulated] = endogenous[is_simulated] + maxe[is_simulated]
 
-    check_prediction_model(endogenous_predicted, model)
-
-    if is_write:
-        record_solution_prediction(results, file_sim)
+    # Checks
+    check_prediction_model(endogenous_predicted, beta)
 
     return predictions
 
 
-def check_prediction_model(predictions_diff, model):
-    """ Perform some basic consistency checks for the prediction model.
-    """
-    # Construct auxiliary object
-    results = model.fit()
-    # Perform basic checks
+def check_prediction_model(predictions_diff, beta):
+    """ Perform some basic consistency checks for the prediction model."""
     assert np.all(predictions_diff >= 0.00)
-    assert results.params.shape == (9,)
-    assert np.all(np.isfinite(results.params))
+    assert beta.shape == (9,)
+    assert np.all(np.isfinite(beta))
 
 
 def calculate_wages_systematic(states, covariates, coeffs_a, coeffs_b, type_shifts):
diff --git a/respy/scripts/scripts_compare.py b/respy/scripts/scripts_compare.py
index ccda415f..9cf4678e 100755
--- a/respy/scripts/scripts_compare.py
+++ b/respy/scripts/scripts_compare.py
@@ -4,12 +4,12 @@ import os
 import shutil
 
 import numpy as np
-from statsmodels.tools.eval_measures import rmse
 
 from respy import RespyCls
 from respy.custom_exceptions import UserError
 from respy.pre_processing.data_processing import process_dataset
 from respy.python.shared.shared_auxiliary import dist_class_attributes
+from respy.python.shared.shared_auxiliary import rmse
 from respy.python.simulate.simulate_auxiliary import construct_transition_matrix
 from respy.python.simulate.simulate_auxiliary import format_float
 from respy.scripts.scripts_update import scripts_update
diff --git a/setup.py b/setup.py
index 55d4e379..dc0cffce 100644
--- a/setup.py
+++ b/setup.py
@@ -91,7 +91,6 @@ def setup_package():
             "numba>=0.43",
             "pandas>=0.24",
             "scipy>=0.19",
-            "statsmodels>=0.9",
             "pytest>=4.0",
             "pyaml",
         ],
diff --git a/tox.ini b/tox.ini
index 7928884b..a6415676 100644
--- a/tox.ini
+++ b/tox.ini
@@ -17,7 +17,6 @@ conda_deps =
     numba >= 0.43
     pandas >= 0.24
     scipy >= 1.2.1
-    statsmodels >= 0.9.0
 conda_channels =
     numba
     defaults

Drop statsmodels dependency
Currently we use statsmodels to calculate the OLS regressions in the value function interpolation and to calculate an rmse in several places. I propose we drop statsmodels for the following reasons:

- Statsmodels causes many deprecation warnings, all of them related to stuff we don't even use
- The least squares problem is solved inefficiently in statsmodels (at least last time I checked). Here are some sources that suggest improvements: [1](https://stackoverflow.com/questions/41648246/efficient-computation-of-the-least-squares-algorithm-in-numpy), [2](https://www.johndcook.com/blog/2010/01/19/dont-invert-that-matrix/) and [3](https://gist.github.com/aldro61/5889795)
- Calculating OLS parameters cannot be done without also calculating standard errors and I don't think we need that for our purpose, so it just adds computation time.
diff --git a/respy/tests/test_f2py.py b/respy/tests/test_f2py.py
index 6e1e8e2a..7aef4a7b 100644
--- a/respy/tests/test_f2py.py
+++ b/respy/tests/test_f2py.py
@@ -5,7 +5,6 @@ from functools import partial
 import numpy as np
 import pytest
 import scipy
-import statsmodels.api as sm
 from numpy.testing import assert_array_almost_equal
 from numpy.testing import assert_array_equal
 from numpy.testing import assert_equal
@@ -24,6 +23,7 @@ from respy.python.shared.shared_auxiliary import extract_cholesky
 from respy.python.shared.shared_auxiliary import get_conditional_probabilities
 from respy.python.shared.shared_auxiliary import get_emaxs_of_subsequent_period
 from respy.python.shared.shared_auxiliary import get_optim_paras
+from respy.python.shared.shared_auxiliary import ols
 from respy.python.shared.shared_auxiliary import read_draws
 from respy.python.shared.shared_auxiliary import replace_missing_values
 from respy.python.shared.shared_auxiliary import transform_disturbances
@@ -212,29 +212,20 @@ class TestClass(object):
             endog = np.dot(exog, beta) + tiny
 
             # Run statsmodels
-            results = sm.OLS(endog, exog).fit()
+            beta_result = ols(y=endog, x=exog)
 
             # Check parameters
-            py = results.params
+            py = beta_result
             f90 = fort_debug.wrapper_get_coefficients(
                 endog, exog, num_covars, num_agents
             )
             assert_almost_equal(py, f90)
 
             # Check prediction
-            py = results.predict(exog)
+            py = exog.dot(beta_result)
             f90 = fort_debug.wrapper_point_predictions(exog, f90, num_agents)
             assert_almost_equal(py, f90)
 
-            # Check coefficient of determination and the standard errors.
-            py = [results.rsquared, results.bse]
-            f90 = fort_debug.wrapper_get_pred_info(
-                endog, f90, exog, num_agents, num_covars
-            )
-
-            for i in range(2):
-                assert_almost_equal(py[i], f90[i])
-
     def test_3(self):
         """ Compare results between FORTRAN and PYTHON of selected functions.
         """
@@ -903,9 +894,7 @@ class TestClass(object):
         )
         assert_almost_equal(endogenous, replace_missing_values(f90))
 
-        py = get_predictions(
-            endogenous, exogenous, max_emax, is_simulated, file_sim, False
-        )
+        py = get_predictions(endogenous, exogenous, max_emax, is_simulated)
 
         f90 = fort_debug.wrapper_get_predictions(
             endogenous,

environment.yml
pip install -e .
parse_log_pytest
environment.yml
pytest pytest-cov pytest-xdist
apt-get update apt-get install -y gfortran libblas-dev liblapack-dev
3.7
pytest --no-header -rA --tb=line --color=no -p no:cacheprovider -W ignore::DeprecationWarning
alabaster==0.7.13
Babel==2.14.0
certifi==2022.12.7
charset-normalizer==3.4.1
codecov==2.1.13
colorama==0.4.6
coverage==7.2.7
distlib==0.3.9
doc8==1.0.0
docutils==0.20.1
exceptiongroup==1.2.2
execnet==2.0.2
filelock==3.12.2
idna==3.10
imagesize==1.4.1
importlib-metadata==6.7.0
iniconfig==2.0.0
jinja2==3.1.6
livereload==2.7.1
llvmlite==0.28.0
MarkupSafe==2.1.5
numba==0.43.1
numpy==1.21.5
packaging==24.0
pandas==0.24.2
patsy==0.5.3
pbr==6.1.1
platformdirs==4.0.0
pluggy==1.2.0
py==1.11.0
pyaml==23.5.8
pygments==2.17.2
pytest==7.4.4
pytest-cov==4.1.0
pytest-xdist==3.5.0
python-dateutil==2.8.2
pytz==2022.7
PyYAML==6.0.1
requests==2.31.0
-e git+https://github.com/OpenSourceEconomics/respy.git@712ac1aa368ba2d4b3aa4928024f6fb40a564212#egg=respy
restructuredtext-lint==1.4.0
ruamel.yaml==0.17.40
ruamel.yaml.clib==0.2.8
scipy==1.2.1
six==1.16.0
snowballstemmer==2.2.0
sphinx==5.3.0
sphinx-autobuild==2021.3.14
sphinx-rtd-theme==2.0.0
sphinxcontrib-applehelp==1.0.2
sphinxcontrib-devhelp==1.0.2
sphinxcontrib-htmlhelp==2.0.0
sphinxcontrib-jquery==4.1
sphinxcontrib-jsmath==1.0.1
sphinxcontrib-qthelp==1.0.3
sphinxcontrib-serializinghtml==1.1.5
statsmodels==0.9.0
stevedore==3.5.2
tomli==2.0.1
tornado==6.2
tox==3.28.0
tox-conda==0.10.2
typing-extensions==4.7.1
urllib3==2.0.7
virtualenv==20.26.6
zipp==3.15.0

name: respy
channels:
  - defaults
  - https://repo.anaconda.com/pkgs/main
  - https://repo.anaconda.com/pkgs/r
  - conda-forge
dependencies:
  - _libgcc_mutex=0.1=main
  - _openmp_mutex=5.1=1_gnu
  - blas=1.0=openblas
  - ca-certificates=2025.2.25=h06a4308_0
  - certifi=2022.12.7=py37h06a4308_0
  - ld_impl_linux-64=2.40=h12ee557_0
  - libffi=3.4.4=h6a678d5_1
  - libgcc-ng=11.2.0=h1234567_1
  - libgfortran-ng=7.5.0=ha8ba4b0_17
  - libgfortran4=7.5.0=ha8ba4b0_17
  - libgomp=11.2.0=h1234567_1
  - libopenblas=0.3.18=hf726d26_0
  - libstdcxx-ng=11.2.0=h1234567_1
  - llvmlite=0.28.0=py37hd408876_0
  - ncurses=6.4=h6a678d5_0
  - numba=0.43.1=py37h962f231_0
  - numpy=1.21.5=py37h7a5d4dd_2
  - numpy-base=1.21.5=py37hb8be1f0_2
  - openssl=1.1.1w=h7f8727e_0
  - pandas=0.24.2=py37he6710b0_0
  - patsy=0.5.3=py37h06a4308_0
  - pip=19.0.3=py37_0
  - python=3.7.16=h7a1cb2a_0
  - python-dateutil=2.8.2=pyhd3eb1b0_0
  - pytz=2022.7=py37h06a4308_0
  - readline=8.2=h5eee18b_0
  - scipy=1.2.1=py37he2b7bc3_0
  - setuptools=65.6.3=py37h06a4308_0
  - six=1.16.0=pyhd3eb1b0_1
  - sqlite=3.45.3=h5eee18b_0
  - statsmodels=0.9.0=py37h035aef0_0
  - tbb=2022.0.0=hdb19cb5_0
  - tk=8.6.14=h39e8969_0
  - wheel=0.38.4=py37h06a4308_0
  - xz=5.6.4=h5eee18b_1
  - zlib=1.2.13=h5eee18b_1
  - pip:
      - alabaster==0.7.13
      - babel==2.14.0
      - charset-normalizer==3.4.1
      - codecov==2.1.13
      - colorama==0.4.6
      - coverage==7.2.7
      - distlib==0.3.9
      - doc8==1.0.0
      - docutils==0.20.1
      - exceptiongroup==1.2.2
      - execnet==2.0.2
      - filelock==3.12.2
      - idna==3.10
      - imagesize==1.4.1
      - importlib-metadata==6.7.0
      - iniconfig==2.0.0
      - jinja2==3.1.6
      - livereload==2.7.1
      - markupsafe==2.1.5
      - packaging==24.0
      - pbr==6.1.1
      - platformdirs==4.0.0
      - pluggy==1.2.0
      - py==1.11.0
      - pyaml==23.5.8
      - pygments==2.17.2
      - pytest==7.4.4
      - pytest-cov==4.1.0
      - pytest-xdist==3.5.0
      - pyyaml==6.0.1
      - requests==2.31.0
      - restructuredtext-lint==1.4.0
      - ruamel-yaml==0.17.40
      - ruamel-yaml-clib==0.2.8
      - snowballstemmer==2.2.0
      - sphinx==5.3.0
      - sphinx-autobuild==2021.3.14
      - sphinx-rtd-theme==2.0.0
      - sphinxcontrib-applehelp==1.0.2
      - sphinxcontrib-devhelp==1.0.2
      - sphinxcontrib-htmlhelp==2.0.0
      - sphinxcontrib-jquery==4.1
      - sphinxcontrib-jsmath==1.0.1
      - sphinxcontrib-qthelp==1.0.3
      - sphinxcontrib-serializinghtml==1.1.5
      - stevedore==3.5.2
      - tomli==2.0.1
      - tornado==6.2
      - tox==3.28.0
      - tox-conda==0.10.2
      - typing-extensions==4.7.1
      - urllib3==2.0.7
      - virtualenv==20.26.6
      - zipp==3.15.0
prefix: /opt/conda/envs/respy

respy/tests/test_f2py.py::TestClass::test_8 respy/tests/test_f2py.py::TestClass::test_12 respy/tests/test_f2py.py::TestClass::test_11 respy/tests/test_f2py.py::TestClass::test_9
respy/tests/test_f2py.py::TestClass::test_3 respy/tests/test_f2py.py::TestClass::test_1 respy/tests/test_f2py.py::TestClass::test_2 respy/tests/test_f2py.py::TestClass::test_7 respy/tests/test_f2py.py::TestClass::test_10 respy/tests/test_f2py.py::TestClass::test_6 respy/tests/test_f2py.py::TestClass::test_4 respy/tests/test_f2py.py::TestClass::test_5
setup.py development/modules/auxiliary_reliability.py respy/scripts/scripts_compare.py environment.yml tox.ini respy/python/solve/solve_auxiliary.py respy/python/shared/shared_auxiliary.py
setup.py development/modules/auxiliary_reliability.py respy/scripts/scripts_compare.py environment.yml tox.ini respy/python/solve/solve_auxiliary.py respy/python/shared/shared_auxiliary.py