In [1]:

    

import numpy as np
from scipy.optimize import minimize, curve_fit
from scipy.misc import derivative
from IPython.display import display, Math, Latex

def bmatrix(a, pref=''):
    lines = str(a).replace('[', '').replace(']', '').splitlines()
    rv = [r'\begin{bmatrix}']
    rv += ['  ' + ' & '.join(l.split()) + r'\\' for l in lines]
    rv +=  [r'\end{bmatrix}']
    return Math(pref + '\n'.join(rv))

def resources(x):
    return Math('K = {:.4f}, L = {:.4f}'.format(*x))


    

In [2]:

    

outer_coefs = np.array([3.2, 5.8, 14.2, 8.7])
inner_coefs = np.array([4, 8, 0, 3])

p = np.array([20, 25, 15, 15])
I = 1000

bounds = [(coef, None) for coef in inner_coefs]

def partial_derivative(f, var, point=[]):
    args = point[:]

    def wraps(x):
        args[var] = x
        return f(args)

    return derivative(wraps, point[var], dx = 1e-6)

def U(x):
    return sum(outer * np.log(x_i - inner) for outer, x_i, inner in zip(outer_coefs, x, inner_coefs))

def solve(args, objective, constraints, name):
    solution = minimize(
        objective,
        args,
        method='SLSQP',
        bounds=bounds,
        constraints=constraints
    )
    display(bmatrix(solution.x, '{} = '.format(name))) 
    return solution.x


    

In [12]:

    

args = np.array([4.2, 8.7, 0.2, 3.7])
objective = lambda x: -U(x)
constraints = {
    'type': 'ineq',
    'fun': lambda x: I - np.dot(p, x)
}
x = solve(args, objective, constraints, 'x')
display(Math("f(x) = {:.4f}".format(-objective(x))))
L = [partial_derivative(U, i, x) / p[i] for i in range(4)]
display(bmatrix(L, 'L = ')) 
#display(Math("L = {:.4f}".format(partial_derivative(U, 0, x) / p[0])))


    

    






$$x = \begin{bmatrix}
  7.38503176 & 12.90909334 & 20.03276307 & 15.27203902\\
\end{bmatrix}$$






    






$$f(x) = 77.5066$$






    






$$L = \begin{bmatrix}
  0.04726691251732973, & 0.04725923588466685, & 0.04725592077647889, & 0.047261909704351034\\
\end{bmatrix}$$

In [13]:

    

args = np.array([4.7325, 8.2082228, 0.7006161,  3.66873595])
objective = lambda h: np.dot(p, h)
constraints = {
    'type': 'ineq',
    'fun': lambda h: U(h) - U(inner_coefs + 1)
}
h = solve(args, objective, constraints, 'h')
display(Math("f(h) = {:.4f}".format(objective(x))))
L = [p[i] / partial_derivative(U, i, h) for i in range(4)]
display(bmatrix(L, 'L = ')) 
#display(Math("L = {:.4f}".format(p[0]/partial_derivative(U, 0, h))))


    

    






$$h = \begin{bmatrix}
  4.29815315 & 8.43232021 & 1.76407106 & 4.08080261\\
\end{bmatrix}$$






    






$$f(h) = 1000.0001$$






    






$$L = \begin{bmatrix}
  1.863457206027938, & 1.863449172125002, & 1.8634553403607232, & 1.8634527689510587\\
\end{bmatrix}$$

In [5]:

    

K = [49920, 45750, 50550, 505750, 47820, 47900, 51900, 45970, 48030, 48100]
L = [10680, 10310, 10680, 10800, 10040, 10420, 10940, 10710, 9900, 9930]
F = [2860, 2940, 2950, 2880, 2510, 2690, 2990, 2800, 3000, 3070]


    

In [6]:

    

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

p0 = [3.155989, 0.68368306, 0.13993322]
coeffs, _ = curve_fit(cobb_douglas, (K, L) , F, p0)
display(Math("F = {:.2f}  K^{{{:.4f}}} L^{{{:.4f}}}".format(coeffs[0], coeffs[1], coeffs[2])))


    

    






$$F = 1675.79  K^{0.0019} L^{0.0559}$$

In [7]:

    

if abs(coeffs[1] + coeffs[2] - 1) < 1e-3:
    print('Постійний прибуток до масштабу')
elif coeffs[1] + coeffs[2] > 1:
    print('Прибуток збільшується із масштабом')
else:
    print('Прибуток зменшується із масштабом')
    
print('Еластичність заміщення для функції Кобба-Дугласа const = 1')


    

    




Прибуток зменшується із масштабом
Еластичність заміщення для функції Кобба-Дугласа const = 1

In [8]:

    

price = 70
w = [100, 100]

def pi(x):
    return np.dot(w, x) - price * cobb_douglas(x, coeffs[0], coeffs[1], coeffs[2])

bounds1 = [
    (0, None),
    (0, None),
]


    

In [9]:

    

constraint = {
    'type': 'ineq', 
    'fun': lambda x: 10 - (x[0] **2 + x[1] ** 2)**5
}

short_solution = minimize(pi, [1, 1], method='SLSQP', bounds=bounds1, constraints=constraint)
display(resources(short_solution.x))
display(Math('profit = {:.4f}'.format(-pi(short_solution.x))))


    

    






$$K = 0.2190, L = 1.2397$$






    






$$profit = 118234.8969$$

In [11]:

    

def price_func(x):
    return -x / 80 + 8310/83

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

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

def wM(x):
    return (wK(x), wL(x))


def monopoly_pi(x):
    q = cobb_douglas(x, coeffs[0], coeffs[1], coeffs[2])
    mw = wM(x)
    return mw[0] * x[0] + mw[1] * x[1] - price_func(q) * q

monopoly_solution = minimize(monopoly_pi, [1, 1], method='SLSQP', bounds=bounds1, constraints=[])

print("Витрати ресурсів: ")
display(resources(monopoly_solution.x))
print("Ціни ресурсів: ")
display(resources((wK(monopoly_solution.x), wL(monopoly_solution.x))))
print("Ціна: ", price_func(cobb_douglas(monopoly_solution.x, coeffs[0], coeffs[1], coeffs[2])))
print("Обсяг продукції: ", cobb_douglas(monopoly_solution.x, coeffs[0], coeffs[1], coeffs[2]))


    

    




Витрати ресурсів: 






    






$$K = 71.5496, L = 350.3503$$






    




Ціни ресурсів: 






    






$$K = 0.7887, L = 6.7588$$






    




Ціна:  70.82780491139616
Обсяг продукції:  2343.414161305175