In [1]:

    

import numpy as np
from matplotlib import pyplot as plt
from scipy.optimize import minimize, curve_fit
from scipy.misc import derivative


    

In [2]:

    

def U(x): return np.log(x[0] - 3) + 3.4*np.log(x[1]) + 8.2*np.log(x[2] - 11) + 6.7*np.log(x[3] - 4)

p = np.array([1, 2, 3, 4])
I = 1000
bounds = [
    (4,  np.inf),
    (1,  np.inf),
    (12, np.inf),
    (5,  np.inf),
]


    

In [3]:

    

def partial_derivative(f, i, p):
    def wrapper(x):
        a = p
        a[i] = x
        return f(a)
 
    return derivative(wrapper, p[i], dx = 1e-6)


    

In [4]:

    

x0 = [15, 15, 15, 15]

x = minimize(lambda x: -U(x), x0, bounds=bounds, constraints=[{'type': 'ineq', 'fun': lambda x: I - x @ p}]).x
Ls = [partial_derivative(U, i, x) / p[i] for i in range(len(p))]
print('Solution: {}'.format(x))
print('Objective: {:.2f}'.format(U(x)))
print('Lagrange multipliers: {}'.format(Ls))


    

    




Solution: [  52.15025268   83.48315766  145.25903298   86.27658577]
Objective: 88.66
Lagrange multipliers: [0.020345780171737715, 0.02036338742072985, 0.020358652837633901, 0.020358161378908335]

In [5]:

    

x0 = [15, 15, 15, 15]
q = [6, 5, 15, 10]

x = minimize(lambda x: p @ x, x0, bounds=bounds, constraints=[{'type': 'ineq', 'fun': lambda x: U(x) - U(q)}]).x
Ls = [p[i] / partial_derivative(U, i, x) for i in range(len(p))]
print('Solution: {}'.format(x))
print('Objective: {:.2f}'.format(p @ x))
print('Lagrange multipliers: {}'.format(Ls))


    

    




Solution: [  5.34364123   3.98424524  17.40581003   7.92557422]
Objective: 97.23
Lagrange multipliers: [2.3436402309049384, 2.3436730828664629, 2.3435886635414058, 2.3436258005304635]

In [6]:

    

K = [5700, 4740, 4390, 5330, 5200, 4160, 5200, 4690, 5890, 4930]
L = [12245, 13340, 13860, 14400,  14000, 11000, 14145, 13900, 15050,  13060]
F = [140330, 120355, 125000, 137330, 121570, 113100, 133000, 126165, 149000, 120950]


    

In [7]:

    

def cobb_douglas(x, a, b, c):
    return a * x[0]**b * x[1]**c

coeffs, _ = curve_fit(cobb_douglas, (K, L), F)
print('Coeffs: {}'.format(coeffs))


    

    




Coeffs: [ 139.02507691    0.62259873    0.16045785]

In [8]:

    

if np.isclose(coeffs[1] + coeffs[2], 1): print('Profit to scale is const.')
elif coeffs[1] + coeffs[2] > 1:          print('Profit is increasing with the scale.')
else:                                    print('Profit is decreasing with the scale.')

print('Replacement elasticity for Cobb-Douglas function is equals to 1')


    

    




Profit is decreasing with the scale.
Replacement elasticity for Cobb-Douglas function is equals to 1

In [9]:

    

price = 70
w = [100, 100]

def pi(x): return np.dot(w, x) - price * cobb_douglas(x, *coeffs)

bounds = [
    (0, np.inf),
    (0, np.inf),
]
constraint = {
    'type': 'ineq', 
    'fun': lambda x: 100000000 - x[0] - x[1]
}


    

In [10]:

    

short_x = minimize(pi, [1, 1], bounds=bounds, constraints=constraint).x
print('Solution: {}'.format(short_x))
print('Profit: {:.2f}'.format(-pi(short_x)))


    

    




Solution: [ 29328607.78334537  21137714.98190209]
Profit: 1449882747.74

In [11]:

    

long_x = minimize(pi, [1, 1], bounds=bounds).x
print('Solution: {}'.format(long_x))
print('Profit: {:.2f}'.format(-pi(long_x)))


    

    




Solution: [ 96659920.15248498  13702650.4375944 ]
Profit: 1697360921.77

In [12]:

    

def price_func(x): return -x / 80 + 200

def wK(x): return 0.025 * x[0] - 1
def wL(x): return 0.025 * x[1] - 2

def monopoly_pi(x):
    q = cobb_douglas(x, *coeffs)
    return (wK(x), wL(x)) @ x - price_func(q) * q

monopoly_x = minimize(monopoly_pi, [1, 1], bounds=bounds).x

print('Resources usage: {}'.format(monopoly_x))
print('Resources prices: {}'.format([wK(monopoly_x), wL(monopoly_x)]))
print('Price: {}'.format(price_func(cobb_douglas(monopoly_x, *coeffs))))
print('Production volume: {}'.format(cobb_douglas(monopoly_x, *coeffs)))


    

    




Resources usage: [ 196.37637898  116.51225428]
Resources prices: [3.9094094746108627, 0.91280635702792123]
Price: 100.17411635964717
Production volume: 7986.070691228226