<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Probit
$$Y_{i}^{*}=X_{i}\theta+e_{i}$$
$$Y_{i}=\begin{cases}
1, & Y_{i}^{}>0,\
0, & Y_{i}^{}\leq0.
\end{cases}$$
Step1:  Generate data
Step2:  Plot the data and the model
Step3:  Maximize log-likelihood
$$l\left(y|x,\theta\right)=\sum_{i=1}^{n}\log f\left(y_{i}|x_{i},\theta\right)=\sum_{i=1}^{n}\left[y_{i}\log\Phi\left(x_{i}\theta\right)+\left(1-y_{i}\right)\log\left[1-\Phi\left(x_{i}\theta\right)\right]\right]$$
Step4:  Plot objective function, true parameter, and the estimate
Step5:  Solve first order conditions
Step6:  Plot first order condition
Step7:  Plot original data and fitted model
<ASSISTANT_TASK:>
Python Code:
import numpy as np
import matplotlib.pylab as plt
import seaborn as sns
np.set_printoptions(precision=4, suppress=True)
sns.set_context('notebook')
%matplotlib inline
Explanation: Probit
$$Y_{i}^{*}=X_{i}\theta+e_{i}$$
$$Y_{i}=\begin{cases}
1, & Y_{i}^{}>0,\
0, & Y_{i}^{}\leq0.
\end{cases}$$
End of explanation
# True parameter
theta = .5
# Sample size
n = int(1e2)
# Independent variable, N(0,1)
X = np.random.normal(0, 1, n)
# Error term, N(0,1)
e = np.random.normal(0, 1, n)
# Sort data for nice plots
X = np.sort(X)
# Unobservable dependent variable
Ys = X * theta + e
# Generate observable binary variable
Y = np.zeros_like(Ys)
Y[Ys > 0] = 1
Explanation: Generate data
End of explanation
plt.figure(figsize=(16,8))
# Unobservables
plt.subplot(2, 1, 1)
plt.plot(X, X * theta, label='True model')
plt.scatter(X[Ys > 0], Ys[Ys > 0], c='red', label='Unobserved > 0')
plt.scatter(X[Ys < 0], Ys[Ys < 0], c='blue', label='Unobserved < 0')
plt.ylabel(r'$Y^*$')
plt.xlabel(r'$X$')
plt.legend()
# Observables
plt.subplot(2, 1, 2)
plt.scatter(X, Y, c=[], lw=2)
plt.ylabel(r'$Y$')
plt.xlabel(r'$X$')
plt.show()
Explanation: Plot the data and the model
End of explanation
import scipy.optimize as opt
from scipy.stats import norm
# Define objective function
def f(theta, X, Y):
    Q = - np.sum(Y * np.log(1e-3 + norm.cdf(X * theta)) + (1 - Y) * np.log(1e-3 + 1 - norm.cdf(X * theta)))
    return Q
# Run optimization routine
theta_hat = opt.fmin_bfgs(f, 0., args=(X, Y))
print(theta_hat)
Explanation: Maximize log-likelihood
$$l\left(y|x,\theta\right)=\sum_{i=1}^{n}\log f\left(y_{i}|x_{i},\theta\right)=\sum_{i=1}^{n}\left[y_{i}\log\Phi\left(x_{i}\theta\right)+\left(1-y_{i}\right)\log\left[1-\Phi\left(x_{i}\theta\right)\right]\right]$$
End of explanation
# Generate data for objective function plot
th = np.linspace(-3., 3., 1e2)
Q = [f(z, X, Y) for z in th]
# Plot the data
plt.figure(figsize=(8, 4))
plt.plot(th, Q, label='Q')
plt.xlabel(r'$\theta$')
plt.axvline(x=theta_hat, c='red', label='Estimated')
plt.axvline(x=theta, c='black', label='True')
plt.legend()
plt.show()
Explanation: Plot objective function, true parameter, and the estimate
End of explanation
from scipy.optimize import fsolve
# Define the first order condition
def df(theta, X, Y):
    return - np.sum(X * norm.pdf(X * theta) * (Y - norm.cdf(X * theta)))
# Solve FOC
theta_hat = fsolve(df, 0., args=(X, Y))
print(theta_hat)
Explanation: Solve first order conditions
End of explanation
# Generate data for the plot
th = np.linspace(-3., 3., 1e2)
Q = np.array([df(z, X, Y) for z in th])
# Plot the data
plt.figure(figsize=(8, 4))
plt.plot(th, Q, label='Q')
plt.xlabel(r'$\beta$')
plt.axvline(x=theta_hat, c='red', label='Estimated')
plt.axvline(x=theta, c='black', label='True')
plt.axhline(y=0, c='green')
plt.legend()
plt.show()
Explanation: Plot first order condition
End of explanation
plt.figure(figsize=(16, 8))
# Unobservables
plt.subplot(2, 1, 1)
plt.plot(X, X * theta, label='True model')
plt.plot(X, X * theta_hat, label='Fitted model')
plt.scatter(X[Ys > 0], Ys[Ys > 0], c='red', label='Unobserved > 0')
plt.scatter(X[Ys < 0], Ys[Ys < 0], c='blue', label='Unobserved < 0')
plt.ylabel(r'$Y^*$')
plt.xlabel(r'$X$')
plt.legend()
# Observables
plt.subplot(2, 1, 2)
plt.scatter(X, Y, c=[], label='Y')
plt.plot(X, norm.cdf(X * theta), label=r'$\Phi(X\theta)$')
plt.plot(X, norm.cdf(X * theta_hat), label=r'$\Phi(X\hat{\theta})$')
plt.ylabel(r'$Y$')
plt.xlabel(r'$X$')
plt.legend()
plt.show()
Explanation: Plot original data and fitted model
End of explanation
<END_TASK>