% 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 theatre 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 ChoiceTheatreResults.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_t, 44, DATA.NP_forced_t);
DATA.Y = DATA.Y([1:36, 38:44],:); 
DATA.Y = reshape(DATA.Y, 43*DATA.NP_forced_t,1);

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

% Calculating travel cost
DATA.TCtime = DATA.Time_t_forced.*max(DATA.Hinc,1500)./(160*DATA.HHpeop); 
DATA.Cost = (DATA.thtick_forced3+0.4637*DATA.Dist_t_forced + DATA.TCtime/3)/100;
DATA.Cost = reshape(DATA.Cost, 44, DATA.NP_forced_t);
DATA.Cost = DATA.Cost([1:36, 38:44],:); 
DATA.Cost = reshape(DATA.Cost, 43*DATA.NP_forced_t,1);

ASC = diag(ones(42,1));
ASC = [ASC; zeros(1, 42)];

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

EstimOpt.NamesA = {'6. Piętro';'Ateneum';'Baj';'Capitol';'Dramatyczny';'Druga Strefa';'Guliwer';'IMKA';'Kamienica';'Komedia';'Komuna';'Konsekwentny';'Kwadrat';'Lalka'; ...
   'Jana Kiepury';'Młyn';'Montownia';'Muzyczny Roma';'Na Woli';'Nowy';'Ochoty';'OCH-Teatr';'Palladium';'Polonia';'Polski';'Powszechny';'Praga';'Projekt Teatr';'Rampa';'Rozmaitości'; ...
   'Sabat';'Scena Prezentacje';'Scena Współczesna';'Studio';'Studio Buffo';'Studio Teatralne Koło';'Syrena';'Teatr Narodowy';'Teatr Wielki';'Teatr Niemożliwy';'Opera Kameralna'; ...
   'Współczesny'; '-TC/100'};
%% ****************************  Options ****************************


INPUT.Y = DATA.Y;
INPUT.Xa = DATA.X;

EstimOpt.NP = DATA.NP_forced_t; % Number of respondents
EstimOpt.NCT = 1; % Number of choice tasks per person 
EstimOpt.NAlt = 43; % Number of alternatives


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

[INPUT, Results, EstimOpt, OptimOpt] = DataCleanDCE(INPUT,EstimOpt);
B_backup = [4.22765266943891;1.54470879150997;-4.06241357657265;1.14092879978803;1.35884419701239;-1.58176184832956;-3.60394245586826;3.33686028180944;2.07543705344920;3.12094564979511;-2.18948593748507;-3.00941521258277;1.68632364460961;-2.62851052307564;-0.0323948732855743;-4.94367777683962;0.108398880342905;3.92278017580954;0.578715383483628;0.360140602146182;-3.18663227420797;2.33475677074371;1.73943233206887;2.61473905853531;-1.55729286341861;1.19724500982906;-3.41796430592019;2.26718955598470;0.504060943132433;2.41930939038187;5.51030853231193;-0.311213716831613;-3.61613118787334;0.0601766323461996;2.00107069232709;-2.13356149783622;2.42413698876191;1.83852770414719;6.24833519148614;-4.77669945516703;0.584048795361388;-0.798469299426041;10.3905402224945];
Results.MNL = MNL(INPUT,Results, EstimOpt, OptimOpt);

%% Calculating Consumer surplus

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

DATA.TCtime2 = DATA.Time_t_nforced(DATA.Time_t_nforced > 0).*max(DATA.Hinc2,1500)./(160*DATA.HHpeop2); 
DATA.Cost2 = (DATA.thtick_nforced3(DATA.thtick_nforced3>0)+0.4637*DATA.Dist_t_nforced(DATA.Dist_t_nforced>0) + DATA.TCtime2/3)/100;
DATA.Cost2 = reshape(DATA.Cost2, 44, 1699);
DATA.Cost2 = DATA.Cost2([1:36, 38:44],:);
DATA.Cost2 = reshape(DATA.Cost2, 43*1699,1);

ASC = diag(ones(42,1));
ASC = [ASC; zeros(1, 42)];

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('ChoiceTheatreResults.mat')