<jupyter_start><jupyter_text>Short Tutorial on the Use of ruspy As opposed to the [other replication notebook](https://github.com/OpenSourceEconomics/ruspy/blob/master/promotion/replication/replication_iskhakov_et_al_2016.ipynb), here we do not actually replicate a whole paper but just an aspect of Rust (1987) in order to show how the key objects and functions for the estimation and the demand function derivation are used in the ruspy package.  Estimation with the Nested Fixed Point AlgorithmHere, we replicate the results in Table IX of Rust (1987) for the buses of group 4. For different data sets, you can either prepare them from raw data yourself (the raw data can be found [here](https://editorialexpress.com/jrust/research.html)) or visit the [data repository of OpenSourceEconomics](https://github.com/OpenSourceEconomics/zurcher-data) for Rust's (1987) data which allows to prepare the original data the way it is needed for the ruspy package. In order to estimate the parameters of the model we first load the data below.<jupyter_code>data = pd.read_pickle("group_4.pkl")
data.head()<jupyter_output><empty_output><jupyter_text>The data is the first ingredient for the ``estimate`` function of ruspy and the second is the initialization dictionairy (have a look at the documentation for the all possible keys). The estimation procedure, here, assumes the model to have the specifications as determined in the key "model_specifications". In order to use the Nested Fixed Point Algorithm, we specify it in the "approach" key of the "optimizer" dictionairy. We further assign the specific optimization algorithm below and prescribe the use of the analytical gradient of the likelihood function.<jupyter_code>init_dict_nfxp = {
    "model_specifications": {
        "discount_factor": 0.9999,
        "number_states": 90,
        "maint_cost_func": "linear",
        "cost_scale": 1e-3,
    },
    "optimizer": {
        "approach": "NFXP",
        "algorithm": "scipy_L-BFGS-B",
        "gradient": "Yes",
    },
}<jupyter_output><empty_output><jupyter_text>The above two elements are now given to the ``estimate`` function. The function solves the model and gives out two separate dictionairies containing the results of the estimation of the transition and cost parameters, respectively. As expected, the results below reveal that we find the same parameter estimates as Rust (1987) in Table IX.<jupyter_code>result_transitions_nfxp, result_fixp_nfxp = estimate(init_dict_nfxp, data)
result_fixp_nfxp
result_transitions_nfxp<jupyter_output><empty_output><jupyter_text>Estimation with Mathematical Programming with Equilibrium ConstraintsWe now estimate the model for the same data set with the same specifications. Here, we have to specify that we want to use "MPEC" in the "approach" key and specify and optimizer from the NLOPT or IPOPT library. We opt for the LD_SLSQP provided by NLOPT and specify several starting values ("params"), stopping tolerances and bounds for the parameters.<jupyter_code>init_dict_mpec = {
    "model_specifications": {
        "discount_factor": 0.9999,
        "number_states": 90,
        "maint_cost_func": "linear",
        "cost_scale": 1e-3,
    },
    "optimizer": {
        "approach": "MPEC",
        "algorithm": "LD_SLSQP",
        "params": np.concatenate(
            (np.full(90, 0.0), np.array([4.0, 1.0]))
        ),
        "gradient": "Yes",
        "set_ftol_abs": 1e-15,
        "set_xtol_rel": 1e-15,
        "set_xtol_abs": 1e-3,
        "set_lower_bounds": np.concatenate((np.full(90, -np.inf), np.full(2, 0.0))),
        "set_upper_bounds": np.concatenate((np.full(90, 50.0), np.full(2, np.inf))),
    },
}<jupyter_output><empty_output><jupyter_text>We hand the above dictionairy to the ``estimate`` function and present the result dictionairies below. We arrive at the same results again. As MPEC also takes the expected values as parameters, the parameter vector in the cost parameter results also contains the estimated expected values at the beginning.<jupyter_code>result_transitions_mpec, result_fixp_mpec = estimate(init_dict_mpec, data)
result_fixp_mpec
result_transitions_mpec<jupyter_output><empty_output><jupyter_text>Deriving the Implied Demand FunctionFrom the estimated parameters above (or generally from any parameters) we can now derive the implied demand function with ruspy as described in Rust (1987). For this we have to specify a dictionairy that describes the grid of replacement cost for which the expected demand is supposed to be calculated. As the demand calculation involves solving a fixed point one needs to specify the stopping tolerance for contraction iterations. Additionally, one passes for how many months, here 12, and how many buses (37) the implied demand is derived. Additionally, one has to supply the cost and transition parameters and pass it in separately to ``get_demand`` function that coordinates the demand derivation.<jupyter_code>demand_dict = {
    "RC_lower_bound": 4,
    "RC_upper_bound": 13,
    "demand_evaluations": 100,
    "tolerance": 1e-10,
    "num_periods": 12,
    "num_buses": 37,
}
demand_params = np.concatenate((result_transitions_nfxp["x"], result_fixp_nfxp["x"]))<jupyter_output><empty_output><jupyter_text>Lastly, the information for which model specification the demand is supposed to be calculated has to be passed in as a dictionairy. If, as in our case, the parameters come from previously run estimation, we can just pass in the initialization dictionairy used for the estimation. From this the ``get_demand`` function extracts the model information by using the "model_specifications" key. This means that is sufficient to pass in a dictionairy that has this key and under which there is the information as requested for the initialization dictionairy. The "optimizer" key is not needed.Below we now run the demand estimation and give out the results for the first few replacement costs.<jupyter_code>demand = get_demand(init_dict_nfxp, demand_dict, demand_params)
demand.head()<jupyter_output><empty_output><jupyter_text>To get a sense of how the overall demand function looks like, we plot the above results across the whole grid of replacement costs that we specified before.<jupyter_code>plt.plot(demand.index.to_numpy(), demand["demand"].astype(float).to_numpy())
plt.xlabel("Replacement Cost")
plt.ylabel("Expected Annual Demand")<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.07735495945102933, 'fraction_numerical': 0.01372426699937617, 'mean_word_length': 4.926987060998152, 'pattern_counts': {'":': 30, '<': 33, '<?xml version=': 0, '>': 33, 'https://': 3, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '162703', 'license': 'MIT', 'n_tokens_mistral': 1895, 'n_tokens_neox': 1778, 'n_words': 733, 'path': 'promotion/replication/replication.ipynb', 'repo_name': 'MaxBlesch/ruspy', 'size': '6412'}