% 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 museum 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 ChoiceMuseumResults.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_m, 29, DATA.NP_forced_m);
DATA.Y = DATA.Y([1:12, 14:29],:);
DATA.Y = reshape(DATA.Y, 28*DATA.NP_forced_m,1);

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

% Calculating travel cost
DATA.TCtime = DATA.Time_m_forced.*max(DATA.Hinc,1500)./(160*DATA.HHpeop); 
DATA.Cost = (DATA.mustick_forced3+0.4637*DATA.Dist_m_forced + DATA.TCtime/3)/100;
DATA.Cost = reshape(DATA.Cost, 29, DATA.NP_forced_m);
DATA.Cost = DATA.Cost([1:12, 14:29],:);
DATA.Cost = reshape(DATA.Cost, 28*DATA.NP_forced_m,1);
ASC = diag(ones(27,1));
ASC = [ASC; zeros(1, 27)];

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

EstimOpt.NamesA = {'Centrum Nauki Kopernik';'Dom Spotkań z Historią';'Kolejnictwa';'Azji i Pacyfiku';'Fryderyka Chopina';'Historii Polski';'Historii Żydów Polskich';'Historyczne m.st. Warszawy';'Jana Pawła II i Prymasa Wyszyńskiego'; ...
    'Literatury';'Łowiectwa i Jeździectwa';'Marii Sklodowskiej Curie';'Narodowe w Warszawie';'Pałac w Wilanowie';'Powstania Warszawskiego';'Sportu i Turystyki';'MSN';'Teatralne';'Techniki';'Łazienki Królewskie'; ...
    'Wojska Polskiego';'CSW';'Ziemi PAN';'ŻIH';'Państwowe Muzeum Archeologiczne';'Państwowe Muzeum Etnograficzne';'Zachęta'; '-TC/100'};


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


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


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

[INPUT, Results, EstimOpt, OptimOpt] = DataCleanDCE(INPUT,EstimOpt);
B_backup = [1.47397836466442;-3.87344245909521;-2.29681511125823;-4.76564865553049;-1.30802348836028;-5.35656456749099;-0.758211075975957;-2.47048624479657;-4.41449282960702;-6.06145672161487;-3.57375555214655;-5.56061171134965;0.685769783099513;0.0741220978155075;0.410930287650571;-3.94011705603889;-3.45361801589947;-4.63372155347012;-1.19323336167075;-0.619879968008276;-0.872045972474507;-3.55069604325219;-3.95223243161925;-3.34828393604411;-3.22763905976710;-2.03128820622224;-2.01418038192546;10.0409409446871];
Results.MNL = MNL(INPUT,Results, EstimOpt, OptimOpt);

%% Calculating Consumer surplus

DATA.HHpeop2 = reshape(DATA.nrmem(:,ones(29,1))', 29*1699,1);
DATA.Hinc2 = reshape(DATA.moneyhom(:,ones(29,1))', 29*1699,1);
DATA.TCtime2 = DATA.Time_m_nforced(DATA.Time_m_nforced > 0).*max(DATA.Hinc2,1500)./(160*DATA.HHpeop2); 
DATA.Cost2 = (DATA.mustick_nforced3(DATA.Time_m_nforced > 0)+0.4637*DATA.Dist_m_nforced(DATA.Dist_m_nforced>0) + DATA.TCtime2/3)/100;
DATA.Cost2 = reshape(DATA.Cost2, 29, 1699);
DATA.Cost2 = DATA.Cost2([1:12, 14:29],:);
DATA.Cost2 = reshape(DATA.Cost2, 28*1699,1);

ASC = diag(ones(27,1));
ASC = [ASC; zeros(1, 27)];
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('ChoiceMuseumResults.mat')