```python
import numpy as np

#extra imports
import sympy as sm
import scipy as sp
from scipy import linalg
from scipy import optimize
from scipy import interpolate
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
%matplotlib inline
```

# Linear regression

Consider the following **linear equation:**

$$y_i = \beta_0 + \beta_1 x_{1,i} + \beta_2 x_{2,i} + \epsilon_i$$

Assume you have access to data of the **independent variables** ($x_{1,i}$, $x_{2,i}$) and the **dependent variable** ($y_i$) for $N$ individuals, where $i$ indexes individuals. The variable $\epsilon_i$ is a mean-zero **stochastic shock**.

Assume the **data generating process** is given by:


```python
def DGP(N):
    
    # a. independent variables
    x1 = np.random.normal(0,1,size=N)
    x2 = np.random.normal(0,1,size=N)
    
    # b. errors
    eps = np.random.normal(0,1,size=N)
    
    extreme = np.random.uniform(0,1,size=N)
    eps[extreme < 0.05] += np.random.normal(-5,1,size=N)[extreme < 0.05]
    eps[extreme > 0.95] += np.random.normal(5,1,size=N)[extreme > 0.95]
    
    # c. dependent variable
    y = 0.1 + 0.3*x1 + 0.5*x2 + eps
    
    return x1, x2, y
```

**The data you have access to is:**


```python
np.random.seed(2020)
x1,x2,y = DGP(10000)
```

**Question 1:** Estimate the vector of coefficients $\mathbf{\beta} = (\beta_0,\beta_1,\beta_2)$ using **ordinary least squares (OLS)** implemented with **matrix algebra** by

$$ \hat{\mathbf{\beta}} = (\mathbf{X}^{\prime}\mathbf{X})^{-1}\mathbf{X}^{\prime}\mathbf{y} $$

where $\mathbf{X}^{\prime}$ is the transpose of $\mathbf{X}$ and

$$\mathbf{y} = 
\pmatrix{ y_1 \cr y_2 \cr  \vdots \cr y_N 
}
, \quad \mathbf{X} = \pmatrix{
1 & x_{1,1} & x_{2,1} \cr 
1 & x_{1,2} & x_{2,2} \cr 
\vdots & \vdots \cr 
1 & x_{1,N} & x_{2,N} 
}$$


```python
# construct beta0 vector of ones
ones = np.ones(10000)

# construct X matrix
X_trans = np.array([ones,x1, x2])
X = X_trans.T

# set up equation and print answer
beta_hat = linalg.inv(X_trans@X)@X_trans@y
print('the estimate of beta is beta_hat =' + f' ({beta_hat[0]:.2f}, {beta_hat[1]:.2f}, {beta_hat[2]:.2f})')  
```

    the estimate of beta is beta_hat = (0.10, 0.29, 0.50)


**Question 2:** Construct a 3D plot, where the data is plotted as scattered points, and the prediction of the model is given by the plane

$$\hat{y}_i = \hat{\beta}_0 + \hat{\beta}_1 x_{1,i} + \hat{\beta}_2 x_{2,i}$$


```python
# predict y_hat
mesh1,mesh2 = np.meshgrid(x1,x2,indexing='ij')
y_plane = beta_hat[0]+beta_hat[1]*mesh1+beta_hat[2]*mesh2

# plot the predicted plane and the data points (be patient...)
plt.style.use('seaborn-whitegrid')

fig_reg = plt.figure(figsize=(12,8))
ax = fig_reg.add_subplot(1,1,1,projection='3d')
ax.scatter(x1,x2,y,s=5)
ax.plot_surface(mesh1,mesh2,y_plane,cmap=cm.jet)

ax.set_title('actual and predicted values')
ax.set_xlabel('$x_1$')
ax.set_ylabel('$x_2$')
ax.set_zlabel('$y$');
```

**Question 3:** Estimate the vector of coefficients $\mathbf{\beta} = (\beta_0,\beta_1,\beta_2)$ using a **numerical solver** to solve the ordinary least square problem, shown below, directly. Compare your results with the matrix algebra results.

$$ \min_{\mathbf{\beta}} \sum^N_{i=1} (y_i - (\beta_0 + \beta_1 x_{1,i} + \beta_2 x_{2,i}) )^2 $$


```python
# define sum of squares function
def sum_of_squares(beta):
    """
    Arguments: 
        beta: vector of beta coefficients
        
    Returns:
        the sum of squares of the data generating process
    
    """
    return np.sum(np.square(y-(beta[0]+beta[1]*x1+beta[2]*x2)))

```


```python
# initial guess
beta_guess = [1,1,1]

#minimize sum of squares
ols = optimize.minimize(sum_of_squares,x0=beta_guess,method='SLSQP')

#unpack betas
beta_0_ols = ols.x[0]
beta_1_ols = ols.x[1]
beta_2_ols = ols.x[2]

print(f' beta0 = {beta_0_ols:.2f}, beta1 = {beta_1_ols:.2}, beta3 = {beta_2_ols:.2f}')
```

     beta0 = 0.10, beta1 = 0.29, beta3 = 0.50



```python
# compare with matrix algebra 
beta_ols = np.array([beta_0_ols,beta_1_ols,beta_2_ols])
np.allclose(beta_hat,beta_ols)
```




    True



From the comparison, we see that the OLS method makes a (numerically) exact replica of the matrix algebra results.

**Question 4:** Estimate the vector of coefficients $\mathbf{\beta} = (\beta_0,\beta_1,\beta_2)$ using **least absolute deviations (LAD)** using a numerical solver to solve the following problem directly: 

$$  \min_{\beta} \sum^N_{i=1} |y_i - (\beta_0 + \beta_1 x_{1,i} + \beta_2 x_{2,i}) | $$

where $|z|$ is the absolute value of $z$.


```python
# define new loss function
def least_abs_dev(beta):
    """
    Arguments: 
        beta: vector of beta coefficients
        
    Returns:
        the least absolute deviations of the data generating process
    
    """
    return np.sum(np.abs(y-(beta[0]+beta[1]*x1+beta[2]*x2)))

#minimize least absolute deviations
lad = optimize.minimize(least_abs_dev,x0=beta_guess,method='SLSQP')

#unpack betas
beta_0_lad = lad.x[0]
beta_1_lad = lad.x[1]
beta_2_lad = lad.x[2]

print(f' beta0 = {beta_0_lad:.2f}, beta1 = {beta_1_lad:.2}, beta3 = {beta_2_lad:.2f}')
```

     beta0 = 0.09, beta1 = 0.31, beta3 = 0.51



```python
# compare with matrix algebra 
beta_lad = np.array([beta_0_lad,beta_1_lad,beta_2_lad])
np.allclose(beta_hat,beta_lad)
```




    False



From the comparison, we see that the LAD method fails to make a (numerically) exact replica of the matrix algebra results.

**Question 5:** Set $N = 50$. Repeat the estimation using the **OLS** and **LAD** methods $K=5000$ times, drawing a new random sample from the data generating process each time. Compare the estimates from each method using histograms. Which method do you prefer? Explain your choice.


```python
# intialize
N = 50
K = 5000

beta0_ols_vec = np.empty(K)
beta1_ols_vec = np.empty(K)
beta2_ols_vec = np.empty(K)

beta0_lad_vec = np.empty(K)
beta1_lad_vec = np.empty(K)
beta2_lad_vec = np.empty(K)


# run monte carlo
for k in range(K):
    # run DGP
    x1,x2,y = DGP(N)
    
    # find betas
    ols = optimize.minimize(sum_of_squares,x0=beta_guess,method='SLSQP')
    lad = optimize.minimize(least_abs_dev,x0=beta_guess,method='SLSQP')
    
    # store betas
    beta0_ols_vec[k] = ols.x[0]
    beta1_ols_vec[k] = ols.x[1]
    beta2_ols_vec[k] = ols.x[2]
    
    beta0_lad_vec[k] = lad.x[0]
    beta1_lad_vec[k] = lad.x[1]
    beta2_lad_vec[k] = lad.x[2]
```


```python
# plot the simulation results
fig_mc, ax_mc = plt.subplots(ncols=3,figsize=(20,5))

ax_mc[0].hist(beta0_ols_vec,bins=50,alpha=0.5,label='OLS',density=True)
ax_mc[0].hist(beta0_lad_vec,bins=50,alpha=0.5, label='LAD',density=True)
ax_mc[0].axvline(0.1,linestyle='dashed',linewidth=2)
ax_mc[0].set_ylabel('Share of beta estimates in bin')
ax_mc[0].legend(loc='upper left')
ax_mc[0].set_xlabel('Realizations of beta0');

ax_mc[1].hist(beta1_ols_vec,bins=50,alpha=0.5,label='OLS',density=True)
ax_mc[1].hist(beta1_lad_vec,bins=50,alpha=0.5, label='LAD',density=True)
ax_mc[1].axvline(0.3,linestyle='dashed',linewidth=2)
ax_mc[1].legend(loc='upper left')
ax_mc[1].set_xlabel('Realizations of beta1');

ax_mc[2].hist(beta2_ols_vec,bins=50,alpha=0.5,label='OLS',density=True)
ax_mc[2].hist(beta2_lad_vec,bins=50,alpha=0.5, label='LAD',density=True)
ax_mc[2].axvline(0.5,linestyle='dashed',linewidth=2)
ax_mc[2].legend(loc='upper left')
ax_mc[2].set_xlabel('Realizations of beta2');

plt.show()
```

The above figures show the realizations of the MC simulation of the DGP. The dashed line indicates the true value of $\beta_i$. Both methods look unbiased, but it seems the consistency of LAD is better than for OLS. Therefore, we would choose LAD as our preferred method in this case. 


```python
len(beta0_ols_vec)
```




    5000



# Durable purchases

Consider a **household** living in two periods.

In the **second period** it gets utility from **non-durable consumption**, $c$, and **durable consumption**, $d+\chi x$:

$$
\begin{aligned}
v_{2}(m_{2},d)&= \max_{c}\frac{(c^{\alpha}(d+\chi x)^{1-\alpha})^{1-\rho}}{1-\rho}\\
\text{s.t.} \\
x &= m_{2}-c \\
c &\in [0,m_{2}]
\end{aligned}
$$

where 

* $m_2$ is cash-on-hand in the beginning of period 2
* $c$ is non-durable consumption
* $d$ is pre-commited durable consumption
* $x = m_2 - c$ is extra durable consumption
* $\rho > 1$ is the risk aversion coefficient
* $\alpha \in (0,1)$ is the utility weight on non-durable consumption
* $\chi \in (0,1)$ implies that extra durable consumption is *less* valuable than pre-comitted durable consumption
* the second constraint ensures the household *cannot* die in debt

The **value function** $v_2(m_2,d)$ measures the household's value of having $m_2$ at the beginning of period 2 with precomitted durable consumption of $d$. The optimal choice of non-durable consumption is denoted $c^{\ast}(m_2,d)$. The optimal extra durable consumption function is $x^{\ast}(m_2,d) = m_2-c^{\ast}(m_2,d)$.

Define the so-called **end-of-period 1 value function** as:

$$
\begin{aligned}
w(a,d)&\equiv\beta\mathbb{E}_{1}\left[v_2(m_2,d)\right]
\end{aligned}
$$

where 

$$
\begin{aligned}
m_2&= (1+r)a+y \\
y &= \begin{cases}
1-\Delta & \text{with prob. }\frac{1}{3}\\
1 & \text{with prob. }\frac{1}{3}\\
1+\Delta & \text{with prob. }\frac{1}{3}
\end{cases}\\
\end{aligned}
$$

and

* $a$ is assets at the end of period 1
* $\beta > 0$ is the discount factor
* $\mathbb{E}_1$ is the expectation operator conditional on information in period 1
* $y$ is income in period 2
* $\Delta \in (0,1)$ is the level of income risk (mean-preserving)
* $r$ is the return on savings

In the **first period**, the household chooses it's pre-comitted level of durable consumption for the next-period,

$$
\begin{aligned}
v_{1}(m_{1})&=\max_{d} w(a,d)\\&\text{s.t.}&\\
a&= m_{1}-d \\
d&\in [0,m_{1}]\\
\end{aligned}
$$

where $m_1$ is cash-on-hand in period 1. The second constraint ensures the household *cannot* borrow. The **value function** $v_1(m_1)$ measures the household's value of having $m_1$ at the beginning of period 1. The optimal choice of pre-committed durable consumption is denoted $d^{\ast}(m_1)$.

The **parameters** and **grids** for $m_1$, $m_2$ and $d$ should be:


```python
# a. parameters
rho = 2
alpha = 0.8
beta = 0.96
r = 0.04
Delta = 0.25
chi = 0.9 # should be 0.9

# b. grids
m1_vec = np.linspace(1e-8,10,100)
m2_vec = np.linspace(1e-8,10,100)
d_vec = np.linspace(1e-8,5,100)
```

**Question 1:** Find and plot the functions $v_{2}(m_{2},d)$, $c^{\ast}(m_2,d)$, and $x^{\ast}(m_2,d)$. Comment.

First, we define the basic functions:


```python
# basic functions
def x(m2,c):
    x=m2-c
    return x

def v2(c,d,m2,chi,alpha,rho):
    return ((c**(alpha)*(d+chi*x(m2,c))**(1-alpha))**(1-rho))/(1-rho)
```

Then, we define the **solution function for period 2:**


```python
# period 2 solution function
def solve_period_2(rho,alpha,chi):

    # a. grids
    v2_grid = np.empty((100,100))
    c_grid = np.empty((100,100))
    x_grid = np.empty((100,100))

    # b. solve for each m2 in grid
    for i,m2 in enumerate(m2_vec):
        for j,d in enumerate(d_vec):

            # i. objective
            obj = lambda c: -v2(c,d,m2,chi,alpha,rho)

            # ii. initial value (consume third)
            x0 = m2/3

            # iii. optimizer
            result = optimize.minimize_scalar(obj,x0,method='bounded',bounds=[1e-8,m2])

            # iv. save results
            v2_grid[i,j] = -result.fun
            c_grid[i,j] = result.x
            x_grid[i,j] = m2-result.x
        
    return x_grid, v2_grid,c_grid
```

Now we can unpack the solution and plot the functions:


```python
# a. unpack solution and construct grids
x_grid,v2_grid,c_grid = solve_period_2(rho,alpha,chi)
m2_grid,d_grid = np.meshgrid(m2_vec,d_vec,indexing='ij')

# b. initialize figure
fig = plt.figure(figsize=(30,5))


#subplots
# subplot value in period 2
ax = fig.add_subplot(1,3,1,projection='3d')
cs = ax.plot_surface(m2_grid,d_grid,v2_grid,cmap=cm.jet)

# add labels
ax.set_title('utility in period 2')
ax.set_xlabel('$m_2$')
ax.set_ylabel('$d$')
ax.set_zlabel('$v$')

# invert xaxis
ax.invert_xaxis()

# add colorbar
fig.colorbar(cs);


# subplot consumption
ax = fig.add_subplot(1,3,2,projection='3d')
cs = ax.plot_surface(m2_grid,d_grid,c_grid,cmap=cm.jet)

# add labels
ax.set_title('non durable consumption in period 2')
ax.set_xlabel('$m_2$')
ax.set_ylabel('$d$')
ax.set_zlabel('$c$')

# invert xaxis
ax.invert_xaxis()

# add colorbar
fig.colorbar(cs);


# subplot extra durables
ax = fig.add_subplot(1,3,3,projection='3d')
cs = ax.plot_surface(m2_grid,d_grid,x_grid,cmap=cm.jet)

# add labels and title
ax.set_title('extra durable consumption in period 2')
ax.set_xlabel('$m_2$')
ax.set_ylabel('$d$')
ax.set_zlabel('$x$')

# invert xaxis
ax.invert_xaxis()

# add colorbar
fig.colorbar(cs);
```

**Comment:** We see from the plof of $x$ that for extra durable good consumption to be positive, the consumer must have a relatively high amount of cash on hand and/or a relatively low amount of initial durable good consumption. This makes good sense, as the marginal utility from consumption in period 2 ($\alpha=0.8)$ is considerably higher than the marginal utility of durable good consumption ($1-\alpha=0.2$). The kink in the x function is, in a sense, mirrored in the consumption function, as the consumer for a relatively high amount of $d$ spends its full cash endowment on consumption.

The plot of $v_2$ is not that easy to interpret. This is due to a property of the utility function (more precisely the "Cobb-Douglas part"), which tends towards negative infinity for very low values of $d$ and $m_2$. Hence, it is difficult to see that there actually is a change in the dark red area as it appears flat. One workaround is to adjust the starting point of the `np.linspace()` vectors m2_vec and d_vec from 1-e8 to e.g. 0.3, such that values very close to zero are excluded and the readability of the plot improves.


```python
# check if all of m2 is used for consumption for d approximately 2.5
np.allclose(m2_vec[99],c_grid[99,50])
```




    True



**Question 2:** Find and plot the functions $v_{1}(m_{1})$ and $d^{\ast}(m_1)$. Comment.

**Hint:** For interpolation of $v_2(m_2,d)$ consider using `interpolate.RegularGridInterpolator([GRID-VECTOR1,GRID-VECTOR2],VALUE-MATRIX,bounds_error=False,fill_value=None)`.


```python
# basic function
def v1(m1,d,Delta,r,beta,v2_interp):
    a = m1-d
    
    # a. v2 value, if low income
    m2_low = (1+r)*a+(1-Delta)
    v2_low = v2_interp([m2_low,d])[0]
    
    # b. v2 value, if mid income
    m2_mid = (1+r)*a+1
    v2_mid = v2_interp([m2_mid,d])[0]
    
    # b. v2 value, if high income
    m2_high = (1+r)*a+(1+Delta)
    v2_high = v2_interp([m2_high,d])[0]
    
    # c. expected v2 value
    v2_exp = 1/3*v2_low+ 1/3*v2_mid+ 1/3*v2_high
    
    v1 = beta*v2_exp
    
    # d. total value
    return v1
```


```python
# a. define solve function
def solve_period_1(alpha,rho,beta,r,Delta,v1,v2_interp):

    # grids
    v1_vec = np.empty(100)
    d1_vec = np.empty(100)
    
    # solve for each m1 in grid
    for i,m1 in enumerate(m1_vec):
        
        # i. objective
        obj = lambda d: -v1(m1,d,Delta,r,beta,v2_interp)
        
        # ii. initial guess (durable consume half)
        x0 = [m1/2]
        
        # iii. optimize
        result = optimize.minimize_scalar(obj,x0,method='bounded', bounds=(1e-8,m1))
        
        # iv. save
        v1_vec[i] = -result.fun
        d1_vec[i] = result.x
     
    return v1_vec,d1_vec

# b. construct interpolator
v2_interp = interpolate.RegularGridInterpolator((m2_vec,d_vec),v2_grid,
                                                bounds_error=False,fill_value=None)
    
# c. solve period 1
v1_vec,d1_vec = solve_period_1(alpha,rho,beta,r,Delta,v1,v2_interp)


# d. plot
fig = plt.figure(figsize=(15,5))

ax = fig.add_subplot(1,2,1)
ax.plot(m1_vec,v1_vec)
ax.set_title('value function period 1')
ax.set_xlabel('$m_1$')
ax.set_ylabel('$v_1$');

ax = fig.add_subplot(1,2,2)
ax.plot(m1_vec,d1_vec)
ax.set_title('durable consumption')
ax.set_xlabel('$m_1$')
ax.set_ylabel('$d$');
```

**Comment:** We see from the plots that both $v_1$ and $d$ are increasing functions of $m_1$. This makes good sense; the more money you have, the more you can consume (both durably and non-durably), which increases the utility of the agent.

Next, consider an **extension** of the model, where there is also a **period 0**. In this period, the household makes a choice whether to stick with the level of durables it has, $z = 0$, or adjust its stock of durables, $z = 1$. If adjusting, the household loses a part of the value of its durable stock; more specificaly it incurs a proportional loss of $\Lambda \in (0,1)$.

Mathematically, the **household problem in period 0** is:

$$
\begin{aligned}
v_{0}(m_{0},d_{0}) &= \max_{z\in\{0,1\}} \begin{cases}
w(m_{0},d_{0}) & \text{if } z = 0\\
v_1(m_0+(1-\Lambda) d_{0}) & \text{if } z = 1\\
\end{cases}\\
\end{aligned}
$$

The **parameters** and **grids** for $m_0$ and $d_0$ should be:


```python
Lambda = 0.2
m0_vec = np.linspace(1e-8,6,100)
d0_vec = np.linspace(1e-8,3,100)
```

**Question 3:** For which values of $m_0$ and  $d_0$ is the optimal choice not to adjust, i.e. $z = 0$? Show this in a plot. Give an interpretion of your results.


```python
# case 1: z = 0
def solve_period_0_nochange(alpha,rho,beta,r,Delta,v1,v2_interp):

    # a. grids
    v0_vec = np.empty(100)
    d0_vec = np.empty(100)
    
    # b. solve for each m1 in grid
    for i,m0 in enumerate(m0_vec):
        
        # i. objective
        obj = lambda d: -v1(m1,d,Delta,r,beta,v2_interp)
        
        # ii. initial guess (durable consume half)
        x0 = [m1/2]
        
        # iii. optimize
        result = optimize.minimize_scalar(obj,x0,method='bounded', bounds=(1e-8,m1))
        
        # iv. save
        v1_vec[i] = -result.fun
        d1_vec[i] = result.x
     
    return v1_vec,d1_vec

# b. construct interpolator
v2_interp = interpolate.RegularGridInterpolator((m2_vec,d_vec),v2_grid,
                                                bounds_error=False,fill_value=None)
    
# c. solve period 1
v1_vec,d1_vec = solve_period_1(alpha,rho,beta,r,Delta,v1,v2_interp)
```


```python
def w_solve_period_0_nochange(m0,d0,Delta,r,beta,v2_interp):
    
    a = m0-d
    
    # a. v2 value, if low income
    m2_low = (1+r)*a+(1-Delta)
    v2_low = v2_interp([m2_low,d])[0]
    
    # b. v2 value, if mid income
    m2_mid = (1+r)*a+1
    v2_mid = v2_interp([m2_mid,d])[0]
    
    # b. v2 value, if high income
    m2_high = (1+r)*a+(1+Delta)
    v2_high = v2_interp([m2_high,d])[0]
    
    # c. expected v2 value
    v2_exp = 1/3*v2_low+ 1/3*v2_mid+ 1/3*v2_high
    
    v1 = beta*v2_exp
    
    # d. total value
    return v1    
    
#def solve_period_0
```

# Gradient descent

Let $\boldsymbol{x} = \left[\begin{array}{c}
x_1 \\
x_2\\
\end{array}\right]$ be a two-dimensional vector. Consider the following algorithm:

**Algorithm:** `gradient_descent()`

**Goal:** Minimize the function $f(\boldsymbol{x})$.

1. Choose a tolerance $\epsilon>0$, a scale factor $ \Theta > 0$, and a small number $\Delta > 0$
2. Guess on $\boldsymbol{x}_0$ and set $n=1$
3. Compute a numerical approximation of the jacobian for $f$ by

    $$
    \nabla f(\boldsymbol{x}_{n-1}) \approx \frac{1}{\Delta}\left[\begin{array}{c}
    f\left(\boldsymbol{x}_{n-1}+\left[\begin{array}{c}
    \Delta\\
    0
    \end{array}\right]\right)-f(\boldsymbol{x}_{n-1})\\
    f\left(\boldsymbol{x}_{n-1}+\left[\begin{array}{c}
    0\\
    \Delta
    \end{array}\right]\right)-f(\boldsymbol{x}_{n-1})
    \end{array}\right]
    $$

4. Stop if the maximum element in $|\nabla f(\boldsymbol{x}_{n-1})|$ is less than $\epsilon$
5. Set $\theta = \Theta$ 
6. Compute $f^{\theta}_{n} = f(\boldsymbol{x}_{n-1} - \theta \nabla f(\boldsymbol{x}_{n-1}))$
7. If $f^{\theta}_{n} < f(\boldsymbol{x}_{n-1})$ continue to step 9
8. Set $\theta = \frac{\theta}{2}$ and return to step 6     
9. Set $x_{n} = x_{n-1} - \theta \nabla f(\boldsymbol{x}_{n-1})$
10. Set $n = n + 1$ and return to step 3

**Question:** Implement the algorithm above such that the code below can run.

**Optimizer function:**


```python
def gradient_descent(f,x0,epsilon=1e-6,Theta=0.1,Delta=1e-8,max_iter=10000): #step 1
    """
    Algorithm to minimize a two dimensional vector function.
    
    Arguments: 
    f (func): function to be minimized
    x0 (vector): 2D initial guess vector
    epsilon (float): tolerance
    Theta (float): scale factor
    Delta (float): small number to be used in the numerical jacobian
    max_iter (float): maximum number of iterations
    
    Returns:
    x (vector): vector of coordinates that minimize the function within the tolerance
    """
    # step 2
    x = x0 # guess
    n = 1 # iterations
    
    # step 3-10
    while n < max_iter:
        
        # save current x as xprev
        xprev = x
    
        # step 3: compute numerical jacobian
        jac_array = np.array([f([x[0]+Delta,x[1]])-f(x),f([x[0],x[1]+Delta])-f(x)])
        num_jacx = 1/Delta*jac_array
    
        # step 4: check convergence
        if abs(np.max(num_jacx)) < epsilon:
            break
            
        else:
            # step 5:
            theta = Theta 
            
            while True:
                
                # step 6: compute f_theta
                f_theta = f(x-theta*num_jacx)
                
                # step 7: check ineq
                if f_theta < f(x):
                    break

                else:
                    # step 8: update theta
                    theta /= 2

        # step 9: update guess
        x = xprev-theta*num_jacx
              
        # step 10: 
        n += 1
                
    return x,n
```

**Test case:**


```python
def rosen(x):
    return (1.0-x[0])**2+2*(x[1]-x[0]**2)**2

x0 = np.array([1.1,1.1])
try:
    x,it = gradient_descent(rosen,x0)
    print(f'minimum found at ({x[0]:.4f},{x[1]:.4f}) after {it} iterations')
    assert np.allclose(x,[1,1])
except:
    print('not implemented yet')
```

    minimum found at (1.0000,1.0000) after 306 iterations


fde1d0190257f778b9b65ba62daa8d30d2df0bbb
ipynb
Jupyter Notebook
examproject/exam_2020.ipynb
NumEconCopenhagen/projects-2020-kongeraekken-1
07310ce0d1cd6c787f7a6616a044b91dc8189343
MIT
2020-03-10T21:09:46.000Z
2020-03-10T21:09:46.000Z
examproject/exam_2020.ipynb
NumEconCopenhagen/projects-2020-kongeraekken-1
07310ce0d1cd6c787f7a6616a044b91dc8189343
MIT
2020-04-18T08:31:52.000Z
2020-05-15T10:51:48.000Z
examproject/exam_2020.ipynb
NumEconCopenhagen/projects-2020-kongeraekken-1
07310ce0d1cd6c787f7a6616a044b91dc8189343
MIT
2021-03-13T10:51:08.000Z
2021-03-20T13:35:20.000Z
Qwen/Qwen-72B
1. YES
2. YES