Repository: OpenSourceEconomics/respy

Issue:
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.

Discussion:

Patch:
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


Test patch:
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,