# -*- coding: utf-8 -*-
"""
Created on Tue Nov 11 21:38:21 2014
@author: Rolf
Value function iteration with linear interpolation and optimization
Not that it is terribly innovative, but it may be helpful to compare different
approaches or programmes, including Python, R, or GAMS.
The model is a simple deterministic dynamic fisheries problem with no fishing
costs and Gordon-Schaefer growth.
"""
from __future__ import division
import numpy as np
from scipy.optimize import fminbound
from scipy import interp
rho = 0.05 # Discount rate
beta = 1/(1+rho) # Discount factor
gridSize = 500 # Number of gridpoints
maxStock = 100.0 # Maximum stock size, or carrying capacity
minState = 0.0 # Minimum stock size, or lower bound state space
maxState = maxStock # Opper bound state space
growthRate = 0.2 # Intrinsic growth rate
numIter = 10 # Number of iterations
# Define state space
stateSpace = np.linspace(minState,maxState,num=gridSize)
# Vectors needed for the algorithm
valueFunction = np.zeros(gridSize) # Value function
nextValue = np.empty(gridSize) # Temporary value function
optH = np.empty(gridSize) # Policy function
lastOptH = np.zeros(gridSize) # Last policy function
convergence = np.zeros(numIter) # Indicator of convergence
# Function to calculate next state from harvest
# Variable tempState will get value later on
nextState = lambda h: max(0.0,tempState+growthRate*tempState*(1-tempState/maxStock)-h)
# Bellman equation
# Note the minus signs; we are using a minimization algorithm to solve a
# maximization problem
def bellMan(h):
if h > tempState or h < 0.0:
out = inf
else:
out = -h-beta*np.interp(nextState(h),stateSpace,valueFunction)
return out
# Iteration
for iter in range(numIter):
for i,tempState in enumerate(stateSpace):
optH[i] = fminbound(bellMan,0.0,tempState-minState) # Find optimal harvest
nextValue[i] = -bellMan(optH[i]) # Calculate value Bellman Equation
convergence[iter] = sum(abs(lastOptH-optH))
lastOptH = np.copy(optH)
valueFunction = np.copy(nextValue) # Update value function
# Optimal steady state stock according to the model
optX1 = np.empty(gridSize)
for i,tempState in enumerate(stateSpace):
optX1[i] = nextState(optH[i])
# Which, for example, would be for the largest 5 stock sizes:
print optX1[495:500]
# Optimal steady state according to the Golden Rule
print (growthRate-rho)/(2*(growthRate/maxStock))