def criterion_ml(theta, data_obs, type_, num_sim):
    """ This function evaluates the criterion for traditional and
    simulation-based maximum likelihood estimation.
    """
    if type_ == 'traditional':
        probs = norm.cdf(theta)
    elif type_ == 'AR-simulator':
        data_sim, _ = simulate_sample(theta, num_sim)
        probs = np.mean(data_sim)
    else:
        raise AssertionError('Invalid request')

    rslt = data_obs * np.log(probs) + (1.0 - data_obs) * np.log(1.0 - probs)

    return np.sum(rslt)

"""
Explanation: Introduction to Structural Estimation
<font size="3"> ... background material available at <a href="https://github.com/softecon/talks">https://github.com/softecon/talks</a>  </font>
Usefulness of Structural Econometric Models


Verification How accurately does the computational model solve the underlying equations of the model for the quantities of interest?


Validation How accurately does the model represent the reality for the quantities of interest?


Uncertainty Quantification How do the various sources of error and uncertainty feed into uncertainty in the model-based prediction of the quantities of interest?


<img src='images/venn_diagram.png' width='750' height='750'>
Goals of Structural Estimation


identify deep parameters and mechanisms governing economic behavior


provide counterfactual policy predictions


determine optimal policies


integrate separate results in unified framework


Deep Parameters


Individual's Utility Function

risk preference
desire for probabilistic protection
elasticity of (intertemporal) substitution
discount factor



Firm's Production Function

elasticity of (input) substitution



Econmic Environment

level of uncertainty



$\rightarrow$ invariant to changes in policy.
Alternative Estimation Strategies


Traditional

Maximum Likelihood  
Method of Moments
Generalized Method of Moments



Simulation-Based

Simulated Maximum Likelihood
Simulated Method of Moments
Indirect Inference



Toy Example
$$
\begin{align}
& \
V^ &= \theta + \epsilon, \qquad \epsilon \sim \mathcal{N}(0, 1)\ 
 & \
D & = \begin{cases}
1 & V^ > 0 \
0 & \text{otherwise}\
\end{cases}
& \
& \
\Pr[D & = 1 \mid \theta]  = P(\theta) = \Phi(\theta)\
\end{align}
$$
<img src='images/distribution_shift.png' width='750' height='750'>
Maximum Likelihood Estimation
We attempt to find the values of the parameters which would have most likely produced the empirical data.
$$
\begin{align}
& \
\hat{\theta}^{ML} = \arg \max_{\theta} \log \mathcal{L}(\theta; D)
\end{align}
$$
In our case, the criterion function is straightforward to derive.
$$
\begin{align}
 \log \mathcal{L}(\theta; D) & = \sum^N_{i} D_i \log P(\theta) + (1 - D_i) \log(1 - P(\theta)) \
\end{align}
$$
Method of Moments
We attempt to minimize the discrepancy between observed and theoretical moments.
$$
\begin{align}
\
\hat{\theta}^{MM} & = \arg \min_{\theta} (\mu(\theta) - m(D))^2 \
& \
\end{align}
$$
For example, we can use the first moments.
$$
\begin{align}
&\ 
\mu(\theta) & = E[D] = P(\theta) = \Phi(\theta) \
m(D)  & = \frac{1}{N} \sum^N_{i = 1} D_i
\end{align}
$$
Generalized Method of Moments
We attempt to minimize a weighted discrepancy between the observed and theoretical moments.
$$
\begin{align}
\
\hat{\theta}^{GMM} = \arg \min_{\theta} (\mu(\theta) - m(D))' W (\mu(\theta) - m(D))
\end{align}
$$


Which and how many moments to choose?


How to choose the weighing matrix?


Simulation-Based Estimation
$$
\begin{align}
\
P(\theta) = \int^{\theta}{-\infty} \phi(\epsilon) d\epsilon = \int^{\infty}{-\infty} I[\theta + \epsilon > 0]\phi(\epsilon) d\epsilon
\end{align}
$$
We could evaluate the choice probability by simply drawing a set of $R$ standard normal deviates and applying the indicator function to each one.  
$$
\begin{align}
\hat{P}(\theta) = \frac{1}{R} \sum^R_{r=1} I[\theta + \epsilon^r > 0]
\end{align}
$$
End of explanation
"""


def criterion_mm(theta, data_obs, type_, num_sim, lambda_):
    """ This function evaluates the criterion for traditional and 
    simulation-based method of moments estimation. 
    """
    if type_ == 'traditional':
        stat = norm.cdf(theta)[0]
    elif type_ == 'AR-simulator':
        data_sim, _ = simulate_sample(theta, num_sim)
        stat = np.mean(data_sim)
    elif type_ == 'smoothed-AR-simulator':
        _, u = simulate_sample(theta, num_sim)
        stat = np.mean(get_smoothed_probabilities(u, lambda_), axis=0)[1]
    else:
        raise AssertionError('Invalid request')

    rslt = (np.mean(data_obs) - stat) ** 2

    return rslt

"""
Explanation: <center> <b><font size='6'> Chatter in Criterion Function </font></b> 
<img src='images/chatter_criterion_function.png' width='750' height='750'></center>
<center> <b><font size='6'> Maximum-Likelihood Estimation, R = 100</font></b> 
<img src='images/criterion_maximum_likelihood_100.png' width='750' height='750'></center>
<center> <b><font size='6'> Method of Moments Estimation, R = 100</font></b> 
<img src='images/criterion_method_moments_100.png' width='750' height='750'></center>
<center> <b><font size='6'> Maximum-Likelihood Estimation, R = 10,000</font></b> 
<img src='images/criterion_maximum_likelihood_10000.png' width='750' height='750'></center>
Logit-Smoothed Accept-Reject Simulator
$$
\begin{align}
\
\hat{P}(\theta) = \frac{1}{R} \sum^R_{r=1} \frac{\exp{(\theta + \epsilon_r) \lambda^{-1}}}{1 + \exp{(\theta + \epsilon_r) \lambda^{-1}}}\
\end{align}
$$


simulated probabilities are smooth in the parameters


no zero-probability events


<center> <b><font size='6'> Logistic Function</font></b> 
<img src='images/logistic_function.png' width='750' height='750'></center>
End of explanation
"""


def criterion_ii(theta, data_obs, num_sim):
    """ This function evaluates the criterion for indirect inference 
    estimation.
    """
    # Run auxiliary model on observed data.
    beta_obs = sm.OLS(data_obs, np.tile(1, len(data_ob))).fit().params[0]

    # Simulate a dataset with candidate parametrization and run auxiliary 
    # model.
    data_sim, _ = simulate_sample(theta, num_sim)
    beta_ii = sm.OLS(data_sim, np.tile(1, len(data_sim))).fit().params[0]

    return (beta_obs - beta_ii) ** 2

"""
Explanation: <center> <b><font size='6'> Method of Moments Estimation, Smoothing</font></b> 
<img src='images/smoothing_functions.png' width='750' height='750'></center>
Indirect Inference
The basic idea is to use an auxiliary model to describe the data and use its parameters as the empirical moments. Our goal then is to minimize the discrepancy between simulated and observed moments. Alternative distance metrics exist and often a weighting matrix is used.
$$
\begin{align}
\
\hat{\theta} = \arg\min_{\theta} (\hat{\beta}^{OBS} -  \hat{\beta}^{II}(\theta))^2
\end{align}
$$
End of explanation
"""


