<jupyter_start><jupyter_text>Hyperbolic discount function: test hypotheses in terms of $\log(k)$This is the second part of our core analyis. Here we quantitatively evaluate the hypotheses about how discounting for the various commodities will change going from fasting to control condition. We will do this by using maximum likeligood parameter estimation for each of the hypotheses, and then compute the $\Delta AIC$ and $\Delta BIC$ metrics to compare the goodness of the models, controling for model complexity (number of parameters). Burnham & Anderson (2002) present heuristics to interpret the $\Delta$ AIC.- $\Delta AIC = 0-2$: litle to distiniguish the models- $\Delta AIC = 4-7$: considerably less support for the model with higher AIC- $\Delta AIC > 10$: essentially no support for the model with higher AIC Kas et al (1995) present similar heuristics to interpret the $\Delta$ BIC.- $\Delta BIC = 0-2$: litle to distiniguish the models- $\Delta BIC = 2-6$: positive evidence against the model with higher BIC- $\Delta BIC = 6-10$: strong evidence against the model with higher BIC- $\Delta BIC > 10$: very strong evidence against the model with higher BIC<jupyter_code>import pandas as pd
import numpy as np
from scipy.stats import cauchy
from scipy.optimize import minimize
def long_to_wide(df, target_param):
    '''Convert long-form data to wide-form data'''
    # pivot long to wide
    df = df.set_index('id')
    df = df.pivot_table(index='id', columns=['commodity', 'condition'], values=target_param)
    # collapse column multi-index 
    df.columns = [' '.join(col).strip() for col in df.columns.values]
    # set index to a column
    df = df.reset_index()
    return df

data_long = pd.read_csv('parameter_estimation_Hyperbolic.csv')
data_wide = long_to_wide(data_long, target_param='logk')

delta_logk_food = (data_wide['food F'] - data_wide['food C']).values
delta_logk_money = (data_wide['money F'] - data_wide['money C']).values
delta_logk_music = (data_wide['music F'] - data_wide['music C']).values<jupyter_output><empty_output><jupyter_text>Create data. Must be in this format as a dict with fields `'delta_food'`, `'delta_money'`, `'delta_music'`.<jupyter_code>data = {'delta_food': delta_logk_food,
        'delta_money': delta_logk_money, 
        'delta_music': delta_logk_music}<jupyter_output><empty_output><jupyter_text>Define our models. We must do it here, where `data` is in the local scope, in order to call the `nll` method as a static method in order to use `scipy.optimize.minimize`.<jupyter_code>n_observations = len(data['delta_food'])


class Model:

    def calc_aic(self):
        return -2*self.ll + 2*self.free_params

    def calc_bic(self):
        return -2*self.ll + np.log(n_observations)*self.free_params

    def fit(self):
        '''Find parameters which minimise the negative log likelihood.'''
        if not self.bounds:
            result = minimize(self.nll, self.x0,
                              method='Nelder-Mead',
                              options={'disp': True})
        else:
            result = minimize(self.nll, self.x0,
                              method='L-BFGS-B', bounds=self.bounds,
                              options={'disp': True})
        self.mlparams = result.x
        self.nll = result.fun
        self.ll = -self.nll
        self.AIC = self.calc_aic()
        self.BIC = self.calc_bic()
        return self


class H1(Model):
    """Our control (trait-only) model which assumes zero change in AUC"""

    name = "1. Trait only"
    x0 = [0.05]
    free_params = len(x0)
    bounds = None

    @staticmethod
    def nll(params):
        return (-sum(cauchy.logpdf(data['delta_food'], loc=0, scale=params[0]) +
                     cauchy.logpdf(data['delta_money'], loc=0, scale=params[0]) +
                     cauchy.logpdf(data['delta_music'], loc=0, scale=params[0])))


class H2(Model):
    """In-domain model"""

    name = "2. In-domain"
    x0 = [+0.25, 0.05]
    free_params = len(x0)
    bounds = None

    @staticmethod
    def nll(params):
        return -sum(cauchy.logpdf(data['delta_food'], loc=params[0], scale=params[1]) +
                   cauchy.logpdf(data['delta_money'], loc=0, scale=params[1]) +
                   cauchy.logpdf(data['delta_music'], loc=0, scale=params[1]))


class H3(Model):
    """Monetary fungibility model"""

    name = "3. Monetary fungibility"
    x0 = [0.25, 0.05]
    free_params = len(x0)
    bounds = None

    @staticmethod
    def nll(params):
        return -sum(cauchy.logpdf(data['delta_food'], loc=params[0], scale=params[1]) +
               cauchy.logpdf(data['delta_money'], loc=params[0], scale=params[1]) +
               cauchy.logpdf(data['delta_music'], loc=0, scale=params[1]))


class H4(Model):
    """Negative spillover model"""

    name = "4. Negative spillover"
    x0 = [0.25, -0.1, 0.05]
    free_params = len(x0)
    bounds = [(0., None), (None, 0.), (0., None)]

    @staticmethod
    def nll(params):
        return -sum(cauchy.logpdf(data['delta_food'], loc=params[0], scale=params[2]) +
               cauchy.logpdf(data['delta_money'], loc=params[1], scale=params[2]) +
               cauchy.logpdf(data['delta_music'], loc=params[1], scale=params[2]))


class H5(Model):
    """Spillover model"""

    name = "5. Spillover"
    x0 = [+0.25, +0.1, 0.05]
    free_params = len(x0)
    bounds = None

    @staticmethod
    def nll(params):
        return -sum(cauchy.logpdf(data['delta_food'], loc=params[0], scale=params[2]) +
               cauchy.logpdf(data['delta_money'], loc=params[1], scale=params[2]) +
               cauchy.logpdf(data['delta_music'], loc=params[1], scale=params[2]))


class H6(Model):
    """State-only model"""

    name = "6. State-only"
    x0 = [+0.25, 0.05]
    free_params = len(x0)
    bounds = None

    @staticmethod
    def nll(params):
        return -sum(cauchy.logpdf(data['delta_food'], loc=params[0], scale=params[1]) +
               cauchy.logpdf(data['delta_money'], loc=params[0], scale=params[1]) +
               cauchy.logpdf(data['delta_music'], loc=params[0], scale=params[1]))<jupyter_output><empty_output><jupyter_text>Create and fit the models.<jupyter_code>models = [H1(), H2(), H3(), H4(), H5(), H6()]

models = [model.fit() for model in models]

for model in models:
    print(f'model: {model.name}, params: {model.mlparams}')<jupyter_output>Optimization terminated successfully.
         Current function value: 365.664764
         Iterations: 21
         Function evaluations: 42
Optimization terminated successfully.
         Current function value: 347.372101
         Iterations: 52
         Function evaluations: 98
Optimization terminated successfully.
         Current function value: 349.159253
         Iterations: 54
         Function evaluations: 105
Optimization terminated successfully.
         Current function value: 340.165061
         Iterations: 125
         Function evaluations: 218
Optimization terminated successfully.
         Current function value: 348.443621
         Iterations: 46
         Function evaluations: 89
model: 1. Trait only, params: [1.0959375]
model: 2. In-domain, params: [1.86257004 1.03288556]
model: 3. Monetary fungibility, params: [0.93814007 1.05183251]
model: 4. Negative spillover, params: [1.86258682 0.         1.03287569]
model: 5. Spillover, params: [1.84629129 0.46764303 0.99219284]
m[...]<jupyter_text>Create and save a useful table comparing the hypotheses. $$W(IC_m) = \frac{-0.5 \Delta IC_m}{\Sigma_i(-0.5 \Delta IC_i)}$$<jupyter_code>def W(x):
    return np.exp(-x/2)

# summarise data in a DataFrame
aic = np.array([model.AIC for model in models])
delta_aic = aic - min(aic)
waic = W(delta_aic)/sum(W(delta_aic))

bic = np.array([model.BIC for model in models])
delta_bic = bic - min(bic)
wbic = W(delta_bic)/sum(W(delta_bic))
info = {'model': [model.name for model in models],
        'n': [model.free_params for model in models],
        'LL': [model.ll for model in models],
        'AIC': [model.AIC for model in models],
        'deltaAIC': delta_aic,
        'wAIC': waic,
        'BIC': [model.BIC for model in models],
        'deltaBIC': delta_bic,
        'wBIC': wbic,
        'parameters': [model.mlparams for model in models]}

results = pd.DataFrame.from_dict(info)
results.to_csv(f'model_comparison_LOGK_hyperbolic.csv', index=False)
results<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1327053773697389, 'fraction_numerical': 0.03600810778585907, 'mean_word_length': 3.196098049024512, 'pattern_counts': {'":': 0, '<': 20, '<?xml version=': 0, '>': 22, 'https://': 0, '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': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '211313', 'license': 'MIT', 'n_tokens_mistral': 2973, 'n_tokens_neox': 2700, 'n_words': 741, 'path': 'analysis/04_analyse_hypotheses_hyperbolic_logk.ipynb', 'repo_name': 'drbenvincent/skrynka-vincent', 'size': '8387'}