<jupyter_start><jupyter_text>Linear regression<jupyter_code>def DGP(N):
    """ data generating process
    
    Args:
        
        N (int): number of observations
        
    Returns:
    
        x1 (ndarray): independent variable x1
        x2 (ndarray): independent variable x2
        y (ndarray): dependent varialbe y
        
    """
    
    # 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]
    
    # d. depenent variable
    y = 0.1 + 0.3*x1 + 0.5*x2 + eps
    
    return x1, x2, y
np.random.seed(2020)
x1,x2,y = DGP(10000)<jupyter_output><empty_output><jupyter_text>Question 1<jupyter_code>def OLS_estimate(x1,x2,y):
    """ compute OLS estimates using matrix algebra
    
    Args:
        
        x1 (ndarray): independent variable x1
        x2 (ndarray): independent variable x2
        y (ndarray): dependent varialbe y
        
    Returns:
    
        betas (ndrarray): estimates
    
        
    """
    
    X = np.vstack((np.ones(x1.size),x1,x2)).T
    betas = (np.linalg.inv(X.T@X)@X.T)@y
    return betas
    
betas = OLS_estimate(x1,x2,y)
for i,beta in enumerate(betas):
    print(f'beta{i} = {beta:.4f}')<jupyter_output>beta0 = 0.0957
beta1 = 0.2929
beta2 = 0.5033<jupyter_text>Question 2<jupyter_code>fig = plt.figure()
ax = fig.add_subplot(1,1,1,projection='3d')

# a. predicted
x1v = np.linspace(x1.min(),x1.max(),100)
x2v = np.linspace(x2.min(),x2.max(),100)
x1g, x2g = np.meshgrid(x1v,x2v,indexing='ij')
yhat = betas[0] + betas[1]*x1g + betas[2]*x2g
ax.plot_wireframe(x1g,x2g,yhat,color='black',alpha=0.1);

# b. scatter
ax.scatter(x1,x2,y,s=50,edgecolor='black',facecolor=colors[0],alpha=0.5);

# c. details
ax.set_xlabel('$x_1$')
ax.set_ylabel('$x_2$')
ax.set_zlabel('$y$')
ax.invert_xaxis()
fig.tight_layout()<jupyter_output><empty_output><jupyter_text>Question 3<jupyter_code>def OLS_objective(betas,x1,x2,y):
    """ OLS objective
    
    Args:
        
        betas (ndarray): current guess on parameters 
        x1 (ndarray): independent variable x1
        x2 (ndarray): independent variable x2
        y (ndarray): dependent varialbe y
        
    Returns:
    
        ressum (ndrarray): sum of squared residuals
        
    """
    
    yhat = betas[0] + betas[1]*x1 + betas[2]*x2
    ressum = np.sum((yhat-y)**2)
    return ressum

def OLS_estimate_alt(x1,x2,y):
    """ compute OLS estimates using numerical optimizer
    
    Args:
        
        x1 (ndarray): independent variable x1
        x2 (ndarray): independent variable x2
        y (ndarray): dependent varialbe y
        
    Returns:
    
        betas (ndarray): parameter estimates 
        
    """    
    
    betas0 = np.array([0.1,0.3,0.5])
    sol = optimize.minimize(OLS_objective,betas0,args=(x1,x2,y),method='Nelder-Mead')
    betas = sol.x
    return betas

betas = OLS_estimate_alt(x1,x2,y)
for i,beta in enumerate(betas):
    print(f'beta{i} = {beta:.4f}')<jupyter_output>beta0 = 0.0957
beta1 = 0.2929
beta2 = 0.5033<jupyter_text>Question 4<jupyter_code>def LAD_objective(betas,x1,x2,y):
    """ LAD objective
    
    Args:
        
        betas (ndarray): current guess on parameters 
        x1 (ndarray): independent variable x1
        x2 (ndarray): independent variable x2
        y (ndarray): dependent varialbe y
        
    Returns:
    
        ressum (ndrarray): sum of absolute residuals
        
    """
    
    yhat = betas[0] + betas[1]*x1 + betas[2]*x2
    ressum = np.sum(np.abs(yhat-y))
    return ressum

def LAD_estimate(x1,x2,y):
    """ compute LAD estimates using numerical optimizer
    
    Args:
        
        x1 (ndarray): independent variable x1
        x2 (ndarray): independent variable x2
        y (ndarray): dependent varialbe y
        
    Returns:
    
        betas (ndarray): parameter estimates
        
    """
    
    betas0 = np.array([0.1,0.3,0.5])
    sol = optimize.minimize(LAD_objective,betas0,args=(x1,x2,y),method='Nelder-Mead')
    betas = sol.x
    return betas
betas = LAD_estimate(x1,x2,y)
for i,beta in enumerate(betas):
    print(f'beta{i} = {beta:.4f}')<jupyter_output>beta0 = 0.0921
beta1 = 0.3074
beta2 = 0.5116<jupyter_text>Question 5<jupyter_code># a. setup
K = 5000
N = 50

# b. allocate
betas_OLS = np.empty((K,3))
betas_LAD = np.empty((K,3))

# c. estimate
for k in range(K):
    x1,x2,y = DGP(N)
    betas_OLS[k,:] = OLS_estimate(x1,x2,y)
    betas_LAD[k,:] = LAD_estimate(x1,x2,y)<jupyter_output><empty_output><jupyter_text>**Mean and std.:**<jupyter_code>for i in range(3):
    print(f'beta{i}')
    print(f' OLS: mean = {np.mean(betas_OLS[:,i]):.4f}, std. = {np.std(betas_OLS[:,i]):.4f}')
    print(f' LAD: mean = {np.mean(betas_LAD[:,i]):.4f}, std. = {np.std(betas_LAD[:,i]):.4f}')<jupyter_output>beta0
 OLS: mean = 0.1026, std. = 0.2749
 LAD: mean = 0.1025, std. = 0.2032
beta1
 OLS: mean = 0.3019, std. = 0.2844
 LAD: mean = 0.3027, std. = 0.2134
beta2
 OLS: mean = 0.5003, std. = 0.2860
 LAD: mean = 0.5008, std. = 0.2116<jupyter_text>**Histograms:**<jupyter_code>fig = plt.figure(figsize=(12,4))

for i in range(3):
    
    ax = fig.add_subplot(1,3,i+1)
    ax.hist(betas_OLS[:,i],bins=100,alpha=0.5,label='OLS');
    ax.hist(betas_LAD[:,i],bins=100,alpha=0.5,label='LAD');
    ax.set_title(f'$\\beta_{{{i}}}$')
    ax.legend(frameon=True)
    
fig.tight_layout()<jupyter_output><empty_output><jupyter_text>Durable consumption<jupyter_code># a. parameters
rho = 2
alpha = 0.8
chi = 0.9
beta = 0.96
r = 0.04
Delta = 0.25

# b. grids
d_vec = np.linspace(1e-8,5,100)
m1_vec = np.linspace(1e-8,10,100)
m2_vec = np.linspace(0.5,10,100)<jupyter_output><empty_output><jupyter_text>Question 1<jupyter_code>def utility(c,d,x,rho,alpha,chi):
    """ utility
    
    Args:
    
        c (float): non-durable consumption
        d (float): pre-commited durable consumption
        x (float): extra durable consumption
        rho (float): CRRA parameter
        alpha (float): utility weight on non-durable consumption
        chi (float): scala factor for extra durable consumption
        
    Returns:
    
        (float): utility of consumption
    
    """
    
    return (c**alpha*(d+chi*x)**(1-alpha))**(1-rho)/(1-rho)

def v2(c,m2,d,rho,alpha,chi):
    """ value of choice in period 2
    
    Args:
    
        c (float): non-durable consumption
        m2 (float): cash-on-hand in beginning of period 2
        d (float): pre-commited durable consumption
        x (float): extra durable consumption
        rho (float): CRRA parameter
        alpha (float): utility weight on non-durable consumption
        chi (float): scala factor for extra durable consumption
        
    Returns:
    
        (float): value-of-choice
    
    """
    
    x = m2-c
    return utility(c,d,x,rho,alpha,chi)
def solve_period_2(m2_vec,d_vec,rho,alpha,chi):
    """ solve consumer problem in period 2 
    
    Args:
    
        m2 (ndarray): vector of cash-on-hand in beginning of period 2
        d (ndarray): vector of pre-commited durable consumption
        d (float): pre-commited durable consumption
        x (float): extra durable consumption
        rho (float): CRRA parameter
        alpha (float): utility weight on non-durable consumption
        chi (float): scala factor for extra durable consumption
        
    Returns:
    
        v2_mat (ndarray): value function in period 2
        c_ast_mat (ndarray): consumption function
        x_ast_mat (ndarray): implied extra durable consumption function
    
    """
    
    # a. allocate
    v2_mat = np.empty((m2_vec.size,d_vec.size)) 
    c_ast_mat = np.empty((m2_vec.size,d_vec.size))
    x_ast_mat = np.empty((m2_vec.size,d_vec.size))
    
    # b. loop over states
    for i,m2 in enumerate(m2_vec):
        for j,d in enumerate(d_vec):
            
            # i. objective
            obj = lambda c: -v2(c,m2,d,rho,alpha,chi)

            # ii. initial value (consume half)
            x0 = m2/2

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

            # iv. save
            v2_mat[i,j] = -result.fun
            c_ast_mat[i,j] = result.x
            x_ast_mat[i,j] = m2-c_ast_mat[i,j]

    return v2_mat,c_ast_mat,x_ast_mat
v2_mat,c_ast_mat,x_ast_mat = solve_period_2(m2_vec,d_vec,rho,alpha,chi)
for ystr,y in [('value',v2_mat),('c',c_ast_mat),('x',x_ast_mat)]:

    fig = plt.figure()
    ax = fig.add_subplot(1,1,1,projection='3d')

    # a. value function
    m2g,dg = np.meshgrid(m2_vec,d_vec,indexing='ij')
    ax.plot_surface(m2g,dg,y,cmap=cm.jet)

    # b. details
    ax.set_title(ystr)
    ax.set_xlabel('$m_1$')
    ax.set_ylabel('$d$')
    ax.set_zlabel('')
    ax.invert_xaxis()
    fig.tight_layout()<jupyter_output><empty_output><jupyter_text>Questions 2<jupyter_code>v2_interp = interpolate.RegularGridInterpolator([m2_vec,d_vec], v2_mat,
                                                bounds_error=False,fill_value=None)   
def w(a,d,r,Delta,v2_interp):
    """ post-decision value function in period 1
    
    Args:
    
        a (float): end-of-period asset
        d (float): pre-commited durable consumption
        r (float): return on savings
        Delta (float): income risk scale factor
        v2_interp (RegularGridInterpolator): interpolator for value function in period 2
        
    Returns:
    
        (ndarray): discounted post-decision value
    
    """
    
    # a. initialize
    w = 0
    y = np.array([1-Delta,1,1+Delta])
    
    # b. loop over shocks
    for i in range(3):
        m2_now = (1+r)*a + y[i]
        v2_now = v2_interp([m2_now,d])[0]
        w += 1/3*v2_now
        
    return beta*w
def v1(d,m1,beta,r,Delta,v2_interp):
    """ post-decision value function in period 1
    
    Args:
    
        d (float): pre-commited durable consumption
        m1 (float): cash-on-hand in the beginning of period 1
        beta (float): discount factor
        r (float): return on savings
        Delta (float): income risk scale factor        
        v2_interp (RegularGridInterpolator): interpolator for value function in period 2
        
    Returns:
    
        (ndarray): value-of-choice
    
    """
    
    a = m1-d
    return w(a,d,r,Delta,v2_interp)
def solve_period_1(m1_vec,beta,r,Delta,v2_interp):
    """ post-decision value function in period 1
    
    Args:
    
        m1 (ndarray): vector cash-on-hand in the beginning of period 1
        beta (float): discount factor
        r (float): return on savings
        Delta (float): income risk scale factor        
        v2_interp (RegularGridInterpolator): interpolator for value function in period 2
        
    Returns:
    
        v1_vec (ndarray): value function in period 1
        d_ast_vec (ndarray): pre-commited durable consumption function
    
    """
    
    # a. grids
    v1_vec = np.empty(m1_vec.size)
    d_ast_vec = np.empty(m1_vec.size)
    
    # b. solve for each m1 in grid
    for i,m1 in enumerate(m1_vec):
        
        # i. objective
        obj = lambda x: -v1(x[0],m1,beta,r,Delta,v2_interp)
        
        # ii. initial guess (pre-commit half)
        x0 = m1*1/3
        
        # iii. optimize
        result = optimize.minimize(obj,[x0],method='L-BFGS-B',bounds=((1e-8,m1),))
        
        # iv. save
        v1_vec[i] = -result.fun
        d_ast_vec[i] = result.x
     
    return v1_vec,d_ast_vec
v1_vec,d_ast_vec = solve_period_1(m1_vec,beta,r,Delta,v2_interp)
fig = plt.figure(figsize=(10,4))
ax = fig.add_subplot(1,2,1)
ax.plot(m1_vec,d_ast_vec)
ax.set_xlabel('$m_1$')
ax.set_ylabel('$d$')
ax.set_title('pre-committed durable consumption')

ax = fig.add_subplot(1,2,2)
ax.plot(m1_vec,v1_vec)
ax.set_xlabel('$m_1$')
ax.set_ylabel('$v_1$')
ax.set_title('value function in period 1')
fig.tight_layout()<jupyter_output><empty_output><jupyter_text>Question 3<jupyter_code>Lambda = 0.2
m0_vec = np.linspace(1e-8,6,100)
d0_vec = np.linspace(1e-8,3,100)
v1_interp = interpolate.RegularGridInterpolator([m1_vec], v1_vec,
                                                bounds_error=False,fill_value=None)   
z = np.empty((m0_vec.size,d_vec.size))
for i,m0 in enumerate(m0_vec):
    for j,d0 in enumerate(d0_vec):
        
        # a. no adjustment
        v_keep = w(m0,d0,r,Delta,v2_interp)
        
        # b. adjustment
        v_adj = v1_interp([m0+(1-Lambda)*d0])[0]
        
        # c. best
        z[i,j] = 0 if v_keep > v_adj else 1
fig = plt.figure()
ax = fig.add_subplot(1,1,1)

m0g,d0g = np.meshgrid(m0_vec,d0_vec,indexing='ij')
ax.scatter(m0g[z == 0],d0g[z == 0],s=4,label='keep')
ax.scatter(m0g[z == 1],d0g[z == 1],s=4,label='adj')

ax.set_xlim([m0_vec[0],m0_vec[-1]])
ax.set_ylim([d0_vec[0],d0_vec[-1]])
ax.set_xlabel('$m_0$')
ax.set_ylabel('$d_0$')
ax.legend(frameon=True)
fig.tight_layout()<jupyter_output><empty_output><jupyter_text>Gradient descent  Question 1<jupyter_code>def gradient_descent(f,x0,epsilon=1e-6,Theta=0.1,Delta=1e-8,max_iter=10_000):
    """ minimize function with gradient descent algorithm
        
    Args:

        f (callable): function
        x0 (np.ndarray): initial guess
        eps (float,optional): tolerance
        Theta (float,optional): initial step-size
        Delta (float,optional): step-size in numerical derivatives
        max_iter (int,optional): maximum number of iterations        
        
    Returns:
    
        x (ndarray): minimum
        it (int): number of iterations used
        
    """
    
    # a. initialize
    x = x0
    fx = f(x0)
    
    # b. iterate
    it = 0
    while it < max_iter:
        
        # i. jacobian
        fp = np.empty(x0.size)
        for i in range(x0.size):
            x_ = x.copy()
            x_[i] = x[i] + Delta                  
            fp[i] = (f(x_)-fx)/Delta
        
        # ii. check convergence
        if np.max(np.abs(fp)) < epsilon: break
        
        # ii. line search
        theta = Theta
        while it < max_iter:
            
            # o. x value
            x_theta = x - theta*fp

            # oo. new function value
            fx_theta = f(x_theta)
            it += 1

            # ooo. break or continue line search
            if fx_theta < fx:
                fx = fx_theta
                x = x_theta
                break
            else:
                theta /= 2
            
    return x,it
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')<jupyter_output>minimum found at (1.0000,1.0000) after 331 iterations
{'fraction_non_alphanumeric': 0.13572145965852025, 'fraction_numerical': 0.04030800133913626, 'mean_word_length': 2.4864612511671336, 'pattern_counts': {'":': 0, '<': 54, '<?xml version=': 0, '>': 51, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '13', 'dir': 'jupyter-structured-clean-dedup', 'id': '298508', 'license': 'MIT', 'n_tokens_mistral': 5946, 'n_tokens_neox': 5680, 'n_words': 1130, 'path': 'exam/solution_2020.ipynb', 'repo_name': 'AskerNC/lectures-2021', 'size': '14935'}