new_contents: # Basic RBC model with full depreciation (Alternate 1)
  #
  # Jesus Fernandez-Villaverde
  # Haverford, July 3, 2013
  #
  # Modified to use a Cython implementation of the inner loop
  import numpy as np
  import math
  import time
  # Import the Cython version of the inner loop
  from innerloop_ext import innerloop
  
  def main_func():
  
      #  1. Calibration
  
      aalpha = 1.0/3.0     # Elasticity of output w.r.t. capital
      bbeta  = 0.95        # Discount factor
  
      # Productivity values
      vProductivity = np.array([0.9792, 0.9896, 1.0000, 1.0106, 1.0212],float)
  
      # Transition matrix
      mTransition   = np.array([[0.9727, 0.0273, 0.0000, 0.0000, 0.0000],
                       [0.0041, 0.9806, 0.0153, 0.0000, 0.0000],
                       [0.0000, 0.0082, 0.9837, 0.0082, 0.0000],
                       [0.0000, 0.0000, 0.0153, 0.9806, 0.0041],
                       [0.0000, 0.0000, 0.0000, 0.0273, 0.9727]],float)
  
      ## 2. Steady State
  
      capitalSteadyState     = (aalpha*bbeta)**(1/(1-aalpha))
      outputSteadyState      = capitalSteadyState**aalpha
      consumptionSteadyState = outputSteadyState-capitalSteadyState
  
      print "Output = ", outputSteadyState, " Capital = ", capitalSteadyState, " Consumption = ", consumptionSteadyState 
  
      # We generate the grid of capital
      vGridCapital           = np.arange(0.5*capitalSteadyState,1.5*capitalSteadyState,0.00001)
  
      nGridCapital           = len(vGridCapital)
      nGridProductivity      = len(vProductivity)
  
      ## 3. Required matrices and vectors
  
      mOutput           = np.zeros((nGridCapital,nGridProductivity),dtype=float)
      mValueFunction    = np.zeros((nGridCapital,nGridProductivity),dtype=float)
      mValueFunctionNew = np.zeros((nGridCapital,nGridProductivity),dtype=float)
      mPolicyFunction   = np.zeros((nGridCapital,nGridProductivity),dtype=float)
      expectedValueFunction = np.zeros((nGridCapital,nGridProductivity),dtype=float)
  
      # 4. We pre-build output for each point in the grid
  
      for nProductivity in range(nGridProductivity):
          mOutput[:,nProductivity] = vProductivity[nProductivity]*(vGridCapital**aalpha)
  
      ## 5. Main iteration
  
      maxDifference = 10.0
      tolerance = 0.0000001
      iteration = 0
  
      log = math.log
      zeros = np.zeros
      dot = np.dot
  
      while(maxDifference > tolerance):
  
          expectedValueFunction = dot(mValueFunction,mTransition.T)
  
          for nProductivity in xrange(nGridProductivity):
  
              # We start from previous choice (monotonicity of policy function)
              gridCapitalNextPeriod = 0
  
              # - Start Inner Loop - #
              mValueFunctionNew, mPolicyFunction = innerloop(bbeta, nGridCapital, gridCapitalNextPeriod, mOutput, nProductivity, vGridCapital, expectedValueFunction, mValueFunction, mValueFunctionNew, mPolicyFunction)
  
              # - End Inner Loop - #
  
          maxDifference = (abs(mValueFunctionNew-mValueFunction)).max()
  
          mValueFunction    = mValueFunctionNew
          mValueFunctionNew = zeros((nGridCapital,nGridProductivity),dtype=float)
  
          iteration += 1
          if(iteration%10 == 0 or iteration == 1):
              print " Iteration = ", iteration, ", Sup Diff = ", maxDifference
  
      return (maxDifference, iteration, mValueFunction, mPolicyFunction)
  
  if __name__ == '__main__':
      # - Start Timer - #
      t1=time.time()
      # - Call Main Function - #
      maxDiff, iterate, mValueF, mPolicyFunction = main_func()
      # - End Timer - #
      t2 = time.time()
      print " Iteration = ", iterate, ", Sup Duff = ", maxDiff
      print " "
      print " My Check = ", mPolicyFunction[1000-1,3-1]
      print " "
      print "Elapse time is ", t2-t1
  
subject: Create RBC_Python_Cython.py