% 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 theatre attendance

% INPUT: 
    % Culture_Valuation.mat - dataset
    % ChoiceTheatreResults.mat - file with the results of choice model
        % estimation (can be generated with Choice_Model_theatre.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 CountTheatreResults.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('ChoiceTheatreResults.mat');
DATA = load('Culture_Valuation.mat');

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

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.hinc/10000).^2, DATA.havech3 == 1,DATA.havech3 == 2, DATA.havech3 == 3, ...
    DATA.Wars/100,(DATA.Wars/100).^2, DATA.havejob,  DATA.edugroup == 2, DATA.edugroup == 3];

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

EstimOpt.NamesA = {'Cons.'; 'Consumer surplus per trip'; 'Household income (10k PLN)' ; 'Household income squared  (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'}; 


%% 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.682331560116526;0.0112016724755532;0.672599118503069;-0.144302891680685;0.162902165023572;0.100649420591593;-0.0488481292235098;-0.770153738378703;1.47089800909889;0.132435703808819;0.196393615683528;0.406327736901998;1.15512891403788;0.495712389156099];
Results.NBP = POISS(INPUT,Results, EstimOpt, OptimOpt);

save('CountTheatreResults.mat')