% 2019-02-05
% "Publicly funded cultural institutions – a comparative economic valuation study"
%     Authors: Aleksandra Wiśniewska, Wiktor Budziński, Mikołaj Czajkowski
%     This code estimates count data model (negative binomial P) on data
%     for museum attendance

% INPUT: 
    % Culture_Valuation.mat - dataset
    % ChoiceMuseumResults.mat - file with the results of choice model
        % estimation (can be generated with Choice_Model_museum.m)
    % Codes uses user-written package available from:
        % https://github.com/czaj/Tools   
    % User-written package for estimation of count data models
        % folder POISS included with these supplementary materials 
    
% OUTPUT
    % Saves CountMuseumResults.mat with the results of estimation
   
clear all
clc

%% Loading and transforming data

global B_backup; % this one is for storing B in case iterations are interrupted with ctrl-c

load('ChoiceMuseumResults.mat');
DATA = load('Culture_Valuation.mat');


DATA.hinc_miss = isnan(DATA.moneyhom) + 0;
DATA.moneyhom(isnan(DATA.moneyhom)) = mean(DATA.moneyhom(~isnan(DATA.moneyhom)));
DATA.hinc = DATA.moneyhom;
DATA.HHpeop = DATA.nrmem;
DATA.Age = 2014 - DATA.yearborn; 
DATA.Wars = 2014 - DATA.yearwars; 

DATA.X = [CS(:,5), DATA.hinc/10000, DATA.havech3 == 1,DATA.havech3 == 2, DATA.havech3 == 3, ...
    DATA.Wars/100,(DATA.Wars/100).^2, DATA.havejob,  DATA.edugroup == 2, DATA.edugroup == 3,  ...
    DATA.HHpeop == 2, DATA.HHpeop == 3, DATA.HHpeop == 4, DATA.HHpeop > 4,DATA.hinc_miss];

EstimOpt.NP = 1699;
EstimOpt.NCT = 1; % Number of choice tasks per person 

INPUT.Xa = [ones(EstimOpt.NP,1), DATA.X];
INPUT.Y = DATA.museum2;

EstimOpt.NamesA = {'Cons.'; 'Consumer surplus per trip'; 'Household income (10k PLN)' ;  ...
     'Have 1 child (no children is base level)'; 'Have 2 children (0 children is base level)'; 'Have 3 children (0 children is base level)'; 'Years living in Warsaw (in 100 years)'; 'Years living in Warsaw squared (in 100 years)';  ...
     'Have job'; 'Middle education'; 'Higher education';  ...
     '2 people hh (base level 1)'; '3 people hh (base level 1)'; '4 people hh (base level 1)'; 'More than 4 people hh (base level 1)';'Household income not reported'};

%% Estimation
INPUT.MissingInd = zeros(1699,1);

[INPUT, Results, EstimOpt, OptimOpt] = DataCleanCDM(INPUT,EstimOpt);

 EstimOpt.NB = 2; % Estimate NBP
 EstimOpt.Censored = 12; % Data are censored at 12
 EstimOpt.RealMin = 1;
 EstimOpt.HessEstFix = 1;

B_backup = [0.507969382625596;0.0321562382053263;0.257584795524294;0.254401296657023;0.00259973925840211;-0.294341496876126;-1.24916654338139;1.94987092018095;0.102091242734992;0.126290505924170;0.304968637260126;-0.0404425687653593;-0.155163876757623;0.0539316601692822;0.101493368263849;-0.109624944503995;1.57252015925684;-0.309473477816143];
Results.NBP = POISS(INPUT,Results, EstimOpt, OptimOpt);

save('CountMuseumResults.mat')