% 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 choice model on data for cinema choice

% INPUT: 
    % Culture_Valuation.mat - dataset
    % ConsumerSurplus_choice.m - function used to calculate welfare measures
    % Codes uses user-written package for estimation of choice models
    % available from:
    %    https://github.com/czaj/DCE
    %    https://github.com/czaj/Tools    
    
% OUTPUT
    % Saves ChoiceCinemaResults.mat with the results of estimation and
    % calculated Consumer surplus

clear all
clc

%% Loading and transforming data

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

DATA = load('Culture_Valuation.mat');

% Transforming choice variable, excluding alternatives which were never
% chosen
DATA.Y = reshape(DATA.ForcedChoice_c, 26, DATA.NP_forced_c);
DATA.Y = DATA.Y([1:12, 14:22, 24:26],:); % Excluding theatre not one attended 
DATA.Y = reshape(DATA.Y, 24*DATA.NP_forced_c,1);

DATA.HHpeop = reshape(DATA.nrmem(DATA.SQ_c == 0,ones(26,1))', 26*sum(DATA.SQ_c == 0),1);
DATA.moneyhom(isnan(DATA.moneyhom)) = mean(DATA.moneyhom(~isnan(DATA.moneyhom)));
DATA.Hinc = reshape(DATA.moneyhom(DATA.SQ_c == 0,ones(26,1))', 26*sum(DATA.SQ_c == 0),1);

% Calculating travel cost
DATA.TCtime = DATA.Time_c_forced.*max(DATA.Hinc,1500)./(160*DATA.HHpeop); 
DATA.Cost = (DATA.cintick_forced3+0.4637*DATA.Dist_c_forced + DATA.TCtime/3)/100;
DATA.Cost = reshape(DATA.Cost, 26, DATA.NP_forced_c);
DATA.Cost = DATA.Cost([1:12, 14:22, 24:26],:); % Excluding theatre not one attended 
DATA.Cost = reshape(DATA.Cost, 24*DATA.NP_forced_c,1);
ASC = diag(ones(23,1));
ASC = [ASC; zeros(1, 23)];

DATA.ASC = repmat(ASC, sum(DATA.SQ_c == 0),1);
DATA.X = [DATA.ASC, -DATA.Cost];

EstimOpt.NamesA = {'Alchemia';'Atlantic';'Cinema City Arkadia';'Cinema City Bemowo';'Cinema City Janki';'Cinema City Mokotów';'Cinema City Promenada';'Cinema City Sadyba';'Femina'; ...
    'Iluzjon';'KC';'Lab';'Kinoteka';'Kultura';'Luna';'Multikinokino Targówek';'Multikinokino Ursynów';'Multikinokino Złote Tarasy';'Multikinokino Wola Park';'Muranów';'Praha'; ...
    'Dom Kultury "Świt"';'Ursus'; '-TC/100'};

%% ****************************  Options ****************************


INPUT.Y = DATA.Y;
INPUT.Xa = DATA.X;
EstimOpt.NP = DATA.NP_forced_c; % Number of respondents
EstimOpt.NCT = 1; % Number of choice tasks per person 
EstimOpt.NAlt = 24; % Number of alternatives


%% **************************** Estimation ****************************

[INPUT, Results, EstimOpt, OptimOpt] = DataCleanDCE(INPUT,EstimOpt);
B_backup = [-4.75641342971482;1.42007770697613;3.33197346054769;2.84385692206640;2.75418821565303;3.26491737688922;3.16267249336845;2.95714018679405;1.09116424205793;-2.43336234922211;-4.85610501415084;-5.37385954694526;1.13434474884725;-2.07588666118759;0.248987845992381;0.954049594938001;4.39678702361960;3.50509900685991;4.25456829154657;-1.16904243520056;0.605569163005527;-4.31273132754729;-1.65202855382748;24.3620096602774];
Results.MNL = MNL(INPUT,Results, EstimOpt, OptimOpt);

%% Calculating Consumer surplus

DATA.HHpeop2 = reshape(DATA.nrmem(:,ones(26,1))', 26*1699,1);
DATA.Hinc2 = reshape(DATA.moneyhom(:,ones(26,1))', 26*1699,1);

DATA.TCtime2 = DATA.Time_c_nforced(DATA.Time_c_nforced > 0).*max(DATA.Hinc2,1500)./(160*DATA.HHpeop2); 
DATA.Cost2 = (DATA.cintick_nforced3(DATA.cintick_nforced3>0)+0.4637*DATA.Dist_c_nforced(DATA.Dist_c_nforced>0) + DATA.TCtime2/3)/100;
DATA.Cost2 = reshape(DATA.Cost2, 26, 1699);
DATA.Cost2 = DATA.Cost2([1:12, 14:22, 24:26],:); 
DATA.Cost2 = reshape(DATA.Cost2, 24*1699,1);

ASC = diag(ones(23,1));
ASC = [ASC; zeros(1, 23)];
DATA.ASC = repmat(ASC, 1699,1);
DATA.X2 = [DATA.ASC, -DATA.Cost2];
EstimOpt.NP = 1699;

[CS] = ConsumerSurplus_choice(DATA.X2, Results.MNL.bhat,Results.MNL.ihess, EstimOpt, 100);

save('ChoiceCinemaResults.mat')