```
from __future__ import division
import numpy as np
import sympy as sym
import matplotlib as mpl
import matplotlib.pyplot as plt
```


```
## Setting the parameters of the model (calibrated values)
R, rho, gamma, eta, L, A, kapa = sym.var('R, rho, gamma, eta, L, A, kapa')
```


```
### Defining variables ###
k, x, y, l, r = sym.var('k, x, y, l, r')

### Defining production function ###
gamma = 0.5
rho = 0.25
f = A*kapa*(eta*(k*x)**rho + (1- eta)*(y*(l/r))**rho)**(gamma/rho) - k*R
```

## Part 1: The first order condition.
Differentiating this function leads us to the first order condition:
$$\frac{(0.5A\eta\kappa(kx)^{1/4}(\eta(kx)^{1/4} - (ly/r)^{1/4}(1 -\eta)))}{k} - R = 0 $$

Solving it directly doesn't work:


```
sym.solve(0.5*A*eta*kapa*(k**(1/4))*(x**(1/4))*(eta*(k*x)**0.25 + (l*y/r)**0.25*(-eta + 1)) - R, k, set=True)
# Double click to enlarge output
```

# Part 2: A closer llok at the polynomial for the FOC

We can rename $K = k^{1/4}$ and obtain a thrid degree polynomial in $K$, and try to solve it by the usual equation solver:


```
## Defining auxiliary functions a, c, d
a = 2*R
c = -A*kapa*(eta**2)*x**0.5
d = -A*kapa*eta*(1-eta)*(x*y*(l/r))**0.25
K = sym.var('K')
pol_K = K**3*a + K*c + d
sym.solve(pol_K, K, warn=True, set=True, rational=True, manual=True)
```




    ([], set())



...which doesn't seem to work. How to solve this cubic equation?

There are several ways. But before proceeding, it is good to visualise the equation for the parameters of the model and reasonable values of the variables.

Next cell plots, for some range of values from 0 to 10, the value of the polynomial we are after. When the line crosses the vertical zero, we have a root. It is easy to see that there are, in fact, 3 real roots, two negatives and one positive.


```
K_range = np.arange(-3.0,3.0,0.01)
polik = sym.lambdify((K,x,y,l,r,A,kapa,eta,R),pol_K)
for i in np.arange(0,10,0.1):
    plt.plot(K_range, polik(K_range,i,i,1,1,1,1,0.89,0.13099))
plt.axvline()
plt.axhline()
plt.show()
```

The graph above shows how for most parameter values, this polynomial has three roots, of which one is positive.

Below is a zoomed in graph for values of $x,y,l$ very close to zero. There are still three roots, only one positive.


```
K_range = np.arange(-1.0,1.0,0.001)
plt.plot(K_range, polik(K_range,0.01,0.01,0.5,1,1,1,0.89,0.13099))
plt.axvline()
plt.axhline()
plt.show()
```

We can also vary $l$ (as we are going to normalise $r$ anyway) and see that it displaces the line to the right or to the left:


```
K_range = np.arange(-2.0,2.0,0.01)
for i in np.arange(0.01,10,0.1):
    plt.plot(K_range, polik(K_range,1,1,i,1,1,1,0.89,0.13099))
plt.axvline()
plt.axhline()
plt.show()
```

The extreme case is when $l=0$, in which case the second root is exactly zero. This is true for any *reasobale* combination of $x$ and $y$:


```
K_range = np.arange(-2.0,2.0,0.01)
for i in np.arange(0.01,10,0.5):
    plt.plot(K_range, polik(K_range,i,i,0,1,1,1,0.89,0.13099))
plt.axvline()
plt.axhline()
plt.show()
```

In the case when $l$ is also close to zero, the thried root gets closer to the zero as well. This makes sense! Zero workers, zero capital.

# Part 3: Solving for K
The first method I was using was to check for the condition below:

if delta < 0 --> only one real root

if delta > 0 --> three real roots

if delta = 0 --> combination of both


```
delta = 4*A*kapa*eta**4*x - 54*R*((1-eta)**2)*((y*(l/r))**(1/2))

def delta_condition(x,y,l,r,A,kappa,eta,R):
    '''
    Condition for one real root
    
    Input: A series of model parameters and values (x,y,l,r,A,kappa,eta,R)
    
    Output: the numerical value of the Delta equation followed to True if Delta<0.
    '''
    delta = 4*A*kappa*(eta**4)*x**(1/2) - 54*R*((1-eta)**2)*((y*(l/r))**(1/2))
    return delta, delta<0

# Example for (1,10,10):
delta_condition(1,1,1,1,1,1,0.89,0.13099)
```




    (2.4241007740000002, False)



As the condition fails, and as we saw in teh grpah above, there are likely to be three real roots.

A way of finding them is the *Trigonometric Method*: If the condition below applies, 


```
def trig_condition(x,y,l,r,A,kappa,eta,R):
    '''
    Condition for three real roots
    
    Input: A series of model parameters and values (x,y,l,r,A,kappa,eta,R)
    
    Output: the numerical value of the trigonometric condition equation followed by the result.
    '''
    arcosine = (27/4)*(1-eta)*(y*l/r)**(1/2) - ((eta**4)*A*kappa*(x**(1/2)))/(2*R)
    if arcosine<=0:
        message = "The equation has three real roots"
    else:
        message = "The equation has at least one imaginary root"
    return arcosine, message
trig_condition(1,1,1,1,1,1,0.89,0.13099)
```




    (-1.6524248415909615, 'The equation has three real roots')



The condition is met! Let's apply the trigonometric solution found on wikipedia for the three roots!


```
p = -(eta**2*A*kapa*(x**0.5))/(2*R)
q = -(eta*(1-eta)*kapa*A*(x*y*(l/r))**(1/4))/(2*R)
qmin = (eta*(1-eta)*kapa*A*(x*y*(l/r))**(1/4))/(2*R)
root_zero = 2*(-p/3)**(0.5)*sym.cos(((1/3)*sym.acos(((3*q)/(2*p))*(-3/p)**(1/2))))
root_two = -(2*(-p/3)**(0.5)*sym.cos(((1/3)*sym.acos(((3*qmin)/(2*p))*(-3/p)**(1/2)))))
root_one =  - root_zero - root_two
print "First root =", root_zero
print "Second root =", root_one
print "Thrird root =", root_two
```

    First root = 2*6**(-0.5)*(A*eta**2*kapa*x**0.5/R)**0.5*cos(0.333333333333333*acos(3*6**0.5*x**(-0.5)*(l*x*y/r)**0.25*(R*x**(-0.5)/(A*eta**2*kapa))**0.5*(-eta + 1)/(2*eta)))
    Second root = 2*6**(-0.5)*(A*eta**2*kapa*x**0.5/R)**0.5*cos(0.333333333333333*acos(-3*6**0.5*x**(-0.5)*(l*x*y/r)**0.25*(R*x**(-0.5)/(A*eta**2*kapa))**0.5*(-eta + 1)/(2*eta))) - 2*6**(-0.5)*(A*eta**2*kapa*x**0.5/R)**0.5*cos(0.333333333333333*acos(3*6**0.5*x**(-0.5)*(l*x*y/r)**0.25*(R*x**(-0.5)/(A*eta**2*kapa))**0.5*(-eta + 1)/(2*eta)))
    Thrird root = -2*6**(-0.5)*(A*eta**2*kapa*x**0.5/R)**0.5*cos(0.333333333333333*acos(-3*6**0.5*x**(-0.5)*(l*x*y/r)**0.25*(R*x**(-0.5)/(A*eta**2*kapa))**0.5*(-eta + 1)/(2*eta)))


Now we neeed to know which one of this is the poisitive root. Let's evaluate it at the values given above.


```
Rt0 = sym.lambdify((x,y,l,r,A,kapa,eta,R), root_zero)
Rt1 = sym.lambdify((x,y,l,r,A,kapa,eta,R), root_one)
Rt2 = sym.lambdify((x,y,l,r,A,kapa,eta,R), root_two)
Rt0(1,1,1,1,1,1,0.89,0.13099), Rt1(1,1,1,1,1,1,0.89,0.13099), Rt2(1,1,1,1,1,1,0.89,0.13099)
```




    (1.7976085400993507, -0.12422961405626354, -1.6733789260430871)



Let's see if we have indeed solve the polynomial:


```
K_range = np.arange(-3.0,3.0,0.01)
plt.plot(K_range, polik(K_range,1,1,1,1,1,1,0.89,0.13099))
plt.plot((Rt0(1,1,1,1,1,1,0.89,0.13099)), 0, 'ro')
plt.plot(Rt1(1,1,1,1,1,1,0.89,0.13099), 0, 'bo')
plt.plot(Rt2(1,1,1,1,1,1,0.89,0.13099), 0, 'go')
plt.axvline()
plt.axhline()
plt.show()
```

Because $k = K^4$, all of the roots will deliver positive capital. However, I have tried all three of them in the solver (see *Examples with k* notebook) and the first (positive) root is the only one that produces strictly positive choices of capital.

That means we have a winner! Below is the expression for the final root:


```
k_star = (2*6**(-0.5)*(A*eta**2*kapa*x**0.5/R)**0.5*sym.cos(0.333333333333333*sym.acos(3*6**0.5*x**(-0.5)*(l*x*y/r)**0.25*(R*x**(-0.5)/(A*eta**2*kapa))**0.5*(-eta + 1)/(2*eta))))**(4)
```

# Appendix: Some comparative statics

Here I adapt the function we found to become plottable:


```
def kfuc2(x,y,l,r,A,kapa,eta,R):
    return (2*6**(-0.5)*(A*eta**2*kapa*x**0.5/R)**0.5*np.cos(0.333333333333333*np.arccos(3*6**0.5*x**(-0.5)*(l*x*y/r)**0.25*(R*x**(-0.5)/(A*eta**2*kapa))**0.5*(-eta + 1)/(2*eta))))**(4)
```

I'm assuming positive assortive matching, that is, high types match with high types.

### Case 1: High types hire more workers

- Labour and capital would be complements:


```
x_range = arange(0.01,1.01,0.01)
plt.plot(x_range, kfuc2(x_range,x_range,x_range*1.8,1,0.5105,1,0.89,0.13099), label="Rich country")
plt.plot(x_range, kfuc2(x_range,x_range,x_range*1.8,1,1,1,0.89,0.83684), label="Poor country")
plt.legend(loc=2)
plt.xlim(0.1,1)
plt.show()
```

### Case 2: High types hire less workers


```
x_range = arange(0.01,1.01,0.01)
plt.plot(x_range, kfuc2(x_range,x_range,(1/x_range)**0.5,1,0.5105,1,0.89,0.13099), label="Rich country")
plt.plot(x_range, kfuc2(x_range,x_range,(1/x_range)**0.5,1,1,1,0.89,0.83684), label="Poor country")
plt.legend(loc=2)
plt.xlim(0.1,1)
plt.show()
```

### Local effects of changing $x,y$ and $l$

- Increasing x,y and l from (1,1,1), rich country


```
x_range = arange(1,1.6,0.01)
plt.plot(x_range, kfuc2(x_range,1,1,1,1,1,0.89,0.13099), label="Increasing $x$")
plt.plot(x_range, kfuc2(1,x_range,1,1,1,1,0.89,0.13099), label="Increasing $y$")
plt.plot(x_range, kfuc2(1,1,x_range,1,1,1,0.89,0.13099), label="Increasing $l$")
plt.legend(loc=2)
plt.xlim(1,1.5)
plt.title("From (1,1,1), the effects on $k^{*}$ of:")
plt.show()
```

- Increasing x,y and l from (1,1,1), poor country (similar change, but less magnitude)


```
x_range = arange(1,1.6,0.01)
plt.plot(x_range, kfuc2(x_range,1,1,1,0.5105,1,0.89,0.83684), label="Increasing $x$")
plt.plot(x_range, kfuc2(1,x_range,1,1,0.5105,1,0.89,0.83684), label="Increasing $y$")
plt.plot(x_range, kfuc2(1,1,x_range,1,0.5105,1,0.89,0.83684), label="Increasing $l$")
plt.legend(loc=2)
plt.xlim(1,1.5)
plt.title("From (1,1,1), the effects on $k^{*}$ of:")
plt.show()
```

- Increasing y and l from 0 to 2, keeping the rest =1 ,rich country


```
x_range = arange(0.01,2,0.01)
plt.plot(x_range, kfuc2(1,x_range,1,1,1,1,0.89,0.13099), label="$y$"+str(i))
plt.plot(x_range, kfuc2(1,1,x_range,1,1,1,0.89,0.13099), label="$l$"+str(i))
plt.show()
```

- Increasing y and l from 0 to 2, keeping the rest =1 ,poor country


```
x_range = arange(0.01,2,0.01)
plt.plot(x_range, kfuc2(1,x_range,1,1,0.5105,1,0.89,0.83684), label="$y$"+str(i))
plt.plot(x_range, kfuc2(1,1,x_range,1,0.5105,1,0.89,0.83684), label="$l$"+str(i))
plt.show()
```

- Increasing $x$ from 0 to 2 and keeping the rest constant at 1 (much bigger, linear effects)


```
plt.plot(x_range, kfuc2(x_range,1,1,1,1,1,0.89,0.13099), label="Rich country")
plt.plot(x_range, kfuc2(x_range,1,1,1,0.5105,1,0.89,0.83684), label="Poor country")
plt.xlim(0.01,2)
plt.legend(loc=2)
plt.show()
```


```

```

12b15a50afa7128b8d321d894dbc3113e84b3bc5
ipynb
Jupyter Notebook
Capital.ipynb
davidrpugh/assortative-matching-large-firms
e475dbc04e59ea066fae681b830fecdb8981c1d6
MIT
2019-07-31T06:34:01.000Z
2020-07-29T10:32:37.000Z
Capital.ipynb
davidrpugh/assortative-matching-large-firms
e475dbc04e59ea066fae681b830fecdb8981c1d6
MIT
Capital.ipynb
davidrpugh/assortative-matching-large-firms
e475dbc04e59ea066fae681b830fecdb8981c1d6
MIT
2016-11-13T19:55:54.000Z
2021-09-17T07:20:22.000Z
Qwen/Qwen-72B
1. YES
2. YES