Merge pull request #74 from davidrpugh/add-solow-model

Add solow model

.gitignore

@@ -8,4 +8,6 @@ build/
 *.TODO
 *.noseids
 *.coverage
-*.h5
+*.h5
+examples/solow_model/depreciation_rates.dta
+examples/solow_model/pwt80.dta

.travis.yml

@@ -25,6 +25,7 @@ install:
   - pip install coveralls coverage
   # - conda install --yes -c dan_blanchard python-coveralls nose-cov
   - python setup.py install
+  - cp quantecon/tests/matplotlibrc .
 
 script:
   - nosetests --with-coverage --cover-package=quantecon

quantecon/models/solow/__init__.py

@@ -0,0 +1,14 @@
+"""
+models directory imports
+
+objects imported here will live in the `quantecon.models.solow` namespace
+
+"""
+__all__ = ['Model', 'CobbDouglasModel', 'CESModel']
+
+from . model import Model
+from . import model
+from . cobb_douglas import CobbDouglasModel
+from . import cobb_douglas
+from . ces import CESModel
+from . import ces

quantecon/models/solow/ces.py

@@ -0,0 +1,124 @@
+"""
+Solow model with constant elasticity of substitution (CES) production.
+
+@author : David R. Pugh
+@date : 2014-12-11
+
+"""
+from __future__ import division
+
+import sympy as sym
+
+from . import model
+
+# declare key variables for the model
+A, k, K, L, Y = sym.symbols('A, k, K, L, Y')
+
+# declare required model parameters
+g, n, s, alpha, delta, sigma = sym.symbols('g, n, s, alpha, delta, sigma')
+
+
+class CESModel(model.Model):
+
+    _required_params = ['g', 'n', 's', 'alpha', 'delta', 'sigma', 'A0', 'L0']
+
+    def __init__(self, params):
+        """
+        Create an instance of the Solow growth model with constant elasticity
+        of subsitution (CES) aggregate production.
+
+        Parameters
+        ----------
+        params : dict
+            Dictionary of model parameters.
+
+        """
+        rho = (sigma - 1) / sigma
+        ces_output = (alpha * K**rho + (1 - alpha) * (A * L)**rho)**(1 / rho)
+        super(CESModel, self).__init__(ces_output, params)
+
+    @property
+    def solow_residual(self):
+        """
+        Symbolic expression for the Solow residual which is used as a measure
+        of technology.
+
+        :getter: Return the symbolic expression.
+        :type: sym.Basic
+
+        """
+        rho = (sigma - 1) / sigma
+        residual = (((1 / (1 - alpha)) * (Y / L)**rho -
+                     (alpha / (1 - alpha)) * (K / L)**rho)**(1 / rho))
+        return residual
+
+    @property
+    def steady_state(self):
+        r"""
+        Steady state value of capital stock (per unit effective labor).
+
+        :getter: Return the current steady state value.
+        :type: float
+
+        Notes
+        -----
+        The steady state value of capital stock (per unit effective labor)
+        with CES production is defined as
+
+        .. math::
+
+            k^* = \left[\frac{1-\alpha}{\bigg(\frac{g+n+\delta}{s}\bigg)^{\rho}-\alpha}\right]^{\frac{1}{rho}}
+
+        where `s` is the savings rate, :math:`g + n + \delta` is the effective
+        depreciation rate, and :math:`\alpha` controls the importance of
+        capital stock relative to effective labor in the production of output.
+        Finally,
+
+        ..math::
+
+            \rho=\frac{\sigma-1}{\sigma}
+
+        where `:math:`sigma` is the elasticity of substitution between capital
+        and effective labor in production.
+
+        """
+        g = self.params['g']
+        n = self.params['n']
+        s = self.params['s']
+        alpha = self.params['alpha']
+        delta = self.params['delta']
+        sigma = self.params['sigma']
+
+        ratio_investment_rates = (g + n + delta) / s
+        rho = (sigma - 1) / sigma
+        k_star = ((1 - alpha) / (ratio_investment_rates**rho - alpha))**(1 / rho)
+
+        return k_star
+
+    def _isdeterminate_steady_state(self, params):
+        """Check that parameters are consistent with determinate steady state."""
+        g = params['g']
+        n = params['n']
+        s = params['s']
+        alpha = params['alpha']
+        delta = params['delta']
+        sigma = params['sigma']
+
+        ratio_investment_rates = (g + n + delta) / s
+        rho = (sigma - 1) / sigma
+
+        return ratio_investment_rates**rho - alpha > 0
+
+    def _validate_params(self, params):
+        """Validate the model parameters."""
+        params = super(CESModel, self)._validate_params(params)
+        if params['alpha'] < 0.0 or params['alpha'] > 1.0:
+            raise AttributeError('Output elasticity must be in (0, 1).')
+        elif params['sigma'] <= 0.0:
+            mesg = 'Elasticity of substitution must be strictly positive.'
+            raise AttributeError(mesg)
+        elif not self._isdeterminate_steady_state(params):
+            mesg = 'Steady state is indeterminate.'
+            raise AttributeError(mesg)
+        else:
+            return params

quantecon/models/solow/cobb_douglas.py

@@ -0,0 +1,106 @@
+"""
+Solow growth model with Cobb-Douglas aggregate production.
+
+@author : David R. Pugh
+@date : 2014-11-27
+
+"""
+from __future__ import division
+
+import numpy as np
+import sympy as sym
+
+from . import model
+
+# declare key variables for the model
+t, X = sym.symbols('t'), sym.DeferredVector('X')
+A, k, K, L = sym.symbols('A, k, K, L')
+
+# declare required model parameters
+g, n, s, alpha, delta = sym.symbols('g, n, s, alpha, delta')
+
+
+class CobbDouglasModel(model.Model):
+
+    _required_params = ['g', 'n', 's', 'alpha', 'delta', 'A0', 'L0']
+
+    def __init__(self, params):
+        """
+        Create an instance of the Solow growth model with Cobb-Douglas
+        aggregate production.
+
+        Parameters
+        ----------
+        params : dict
+            Dictionary of model parameters.
+
+        """
+        cobb_douglas_output = K**alpha * (A * L)**(1 - alpha)
+        super(CobbDouglasModel, self).__init__(cobb_douglas_output, params)
+
+    @property
+    def steady_state(self):
+        r"""
+        Steady state value of capital stock (per unit effective labor).
+
+        :getter: Return the current steady state value.
+        :type: float
+
+        Notes
+        -----
+        The steady state value of capital stock (per unit effective labor)
+        with Cobb-Douglas production is defined as
+
+        .. math::
+
+            k^* = \bigg(\frac{s}{g + n + \delta}\bigg)^\frac{1}{1-\alpha}
+
+        where `s` is the savings rate, :math:`g + n + \delta` is the effective
+        depreciation rate, and :math:`\alpha` is the elasticity of output with
+        respect to capital (i.e., capital's share).
+
+        """
+        s = self.params['s']
+        alpha = self.params['alpha']
+        return (s / self.effective_depreciation_rate)**(1 / (1 - alpha))
+
+    def _validate_params(self, params):
+        """Validate the model parameters."""
+        params = super(CobbDouglasModel, self)._validate_params(params)
+        if params['alpha'] <= 0.0 or params['alpha'] >= 1.0:
+            raise AttributeError('Output elasticity must be in (0, 1).')
+        else:
+            return params
+
+    def analytic_solution(self, t, k0):
+        """
+        Compute the analytic solution for the Solow model with Cobb-Douglas
+        production technology.
+
+        Parameters
+        ----------
+        t : numpy.ndarray (shape=(T,))
+            Array of points at which the solution is desired.
+        k0 : (float)
+            Initial condition for capital stock (per unit of effective labor)
+
+        Returns
+        -------
+        analytic_traj : ndarray (shape=t.size, 2)
+            Array representing the analytic solution trajectory.
+
+        """
+        s = self.params['s']
+        alpha = self.params['alpha']
+
+        # lambda governs the speed of convergence
+        lmbda = self.effective_depreciation_rate * (1 - alpha)
+
+        # analytic solution for Solow model at time t
+        k_t = (((s / (self.effective_depreciation_rate)) * (1 - np.exp(-lmbda * t)) +
+                k0**(1 - alpha) * np.exp(-lmbda * t))**(1 / (1 - alpha)))
+
+        # combine into a (T, 2) array
+        analytic_traj = np.hstack((t[:, np.newaxis], k_t[:, np.newaxis]))
+
+        return analytic_traj