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

% INPUT: 
    % Culture_Valuation.mat - dataset
    % ChoiceCinemaResults.mat - file with the results of choice model
        % estimation (can be generated with Choice_Model_cinema.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 CountCinemaResults.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('ChoiceCinemaResults.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.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,  ...
    DATA.Age/100,DATA.Wars == DATA.Age, DATA.HHpeop == 2, DATA.HHpeop == 3, DATA.HHpeop == 4, DATA.HHpeop > 4];


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

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

%EstimOpt.NamesA = {'Cons.'; 'CS';'HINC/1000' ; 'HINC miss'; 'No. of people';  'No. of kids'; 'Age/100'; 'Years in Wars/100'; 'Female'; 'Have job'; 'Middle education'; 'Higher education'};
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';  ...
      'Age (in 100 years)'; 'Born in Warsaw'; '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)'};


%% 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 = [1.74670423299023;0.0246486226954892;0.731519974469225;-0.180516597817666;0.0352996339777527;-0.140370881152513;-0.238154900321260;-1.45020077797361;2.00373364412791;0.225029655961374;0.186019634981633;0.407293098575221;-1.17522945231850;0.120110524442227;-0.134189491471342;-0.206966801248033;-0.0716731381610988;-0.114023880643333;1.35821023769240;0.682059710190797];
Results.NBP = POISS(INPUT,Results, EstimOpt, OptimOpt);

save('CountCinemaResults.mat')