clear all; close all; clc; 

% Created by Ji Yun Lee (Oct 10th, 2018)
% Revised by Jie Zhao (2019 Mar)

% Paper: journal paper for catastrophe risk insurance
% Case Study: Windstorm Insurance in miami county, FL

% Phase 1: Ex-Ante Individual Homeowner's Decision Model (different with the earthquake case study,
% this code is used to generate the basic information (e.g., wealth, buiding value...) instead of 
% the insurance purchase behaviors, due to the mandatory of this insurance )

% NOTE: All monetary units are thousand dollars. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% INPUT PARAMETERS

% Number of buildings in each type (WSF1, WSF2, WMUH1, WMUH2, MSF1, MSF2,
% CERBL): row (census tract) column (building type) 
TypeBuilding = [264 113 21 6 876 179 76; 252 108 27 11 838 172 73; 235 ...
    101 42 13 781 160 68; 301 129 55 26 999 205 87; 220 94 39 14 732 ...
    150 64; 213 91 31 31 706 145 63; 177 76 44 72 586 120 56];
% Total number of buildings 
Num = sum(TypeBuilding(:)); 
% Total house values: row (census tract) column (building type) 
% Total dollar exposure of each type of building in each census tract 
% Later, it should be divided by the number of buildings in type i in
% census tract j
TotalHouseValue = [70757	30325	13773	2904	234917	48115	20251
70160	30068	21475	7647	232930	47708	20281
73994	31711	48665	24577	245660	50316	21819
74247	31820	66748	23606	246498	50487	21830
58462	25055	53514	21563	194095	39754	17203
45882	19663	22518	24855	152327	31200	14143
38071	16317	46002	115655	126397	25889	15934];
% Building (exterior) values: row (census tract) column (building type) 
% Only exterior house values except content values 
TotalBldgValue = [47172	20217	6732	1492	156613	32077	13501
46773	20045	10648	4515	155286	31805	13521
49330	21141	23272	15234	163775	33544	14546
49498	21213	32658	13320	164332	33658	14554
38975	16704	25164	10112	129399	26503	11469
30588	13109	11330	16068	101552	20800	9429
25381	10878	25110	76174	84265	17259	10623];
% Historical Indices 
Hindex_2018 = 100; Hindex_2010 = 85; 
% House value for each type of building in each census tract: row (census
% tract) column (building type) 
for i = 1:size(TotalHouseValue,1)
    for j = 1:size(TotalHouseValue,2)
        HouseValue_2010(i,j) = TotalHouseValue(i,j)/TypeBuilding(i,j); 
    end
end
% Building value for each type of building in each census tract: row
% (census tract) column (building type)
for i = 1:size(TotalBldgValue,1)
    for j = 1:size(TotalBldgValue,2)
        BldgValue_2010(i,j) = TotalBldgValue(i,j)/TypeBuilding(i,j); 
    end
end
HouseValue = HouseValue_2010*Hindex_2018/Hindex_2010; 
BldgValue = BldgValue_2010*Hindex_2018/Hindex_2010; 
% Return periods 
T = [10 20 50 100 200 500 1000]; 
% Wind speed (mph)
v = [86 104 121 132 142 154 163]; 
% Expected losses: row (type of buildings) column (wind speed, IM) 
% The values are the sum of the total esimated losses of all buildings in
% each building type in all census tracts. It is because this will be only
% used for the purpose of constructing the loss functions. Here, the
% total estimated losses = total - content - inventory (HAZUS-MH) 
Total_EL = [1134	3851	16282	38101	69142	124436	167515
570	2794	16365	35560	58161	90225	108683
135	1684	7996	20426	34431	58776	78293
517	3363	11229	25585	43613	71763	94628
3763	12789	53308	123341	224120	403615	544904
906	4346	25328	55143	90473	141451	170804
62	531	4603	10829	17707	28573	35869];
% Income Range (unit: thousand dollars) 
IncomeRange = [0,10,20,30,40,50,60,75,100,max(max(BldgValue))]; 
% Number of people in the income range: row (census tract) column (income range) 
Income = [50	233	247	163	112	168	143	174	176
97	166	332	64	129	154	160	182	178
68	195	255	84	276	133	182	113	237
227	236	343	186	199	128	179	156	134
148	143	220	173	80	129	90	104	170
83	305	276	214	122	134	164	83	105
490	440	396	404	220	251	31	31	29]; 
% C.O.V. of a thershold of concern (p*) 
COV_p = 0.3; % mean of p* is the same value as the annula probability of occurrence 
% lower and upper bound of a threshold of concern (p*)
lower_p = 1/650; upper_p = 1; 
% Alpha of power utility function 
min_risk = 0.0001; max_risk = 1; mean_risk = 0.5; cov_risk = 0.3; 
% Inflation 
Inf_2018 = 116.5; Inf_2010 = 100; 
% Deductible rate 
D = 0.05; 
% insurance premium per $1,000 of coverage
z = 8.13*10^-3; 
% Minimum coverage 
% The insurance coverage is within the range of [$5,000,the house value]
% from Kunreuther's paper. The inflation-adjusted minimum insurance
% coverage is $6,423 in 2018 dollars. 
min_C = 6.423; % max_C depends on the house value
% Number of intervals of insurance coverage 
N_int_C = 50; 
% Minimum and maximum hazard range 
min_H = 60; max_H = 200; 
% Number of intervals of hazard intensity 
N_int_H = 10; 
% Maximum government assistance
GA_max = 33.792; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ZERO MATRICES 
% Index matrix: row(number of bldgs) column (building type) the first column
% indicates census tract and the second column indicates building type
IndexBldg = zeros(Num,2); 
% A threshold of concern matrix: length(total number of houses)
p_star = zeros(Num,1);
% Alpha of power untiltiy function matrix: length(total nubmer of houses)
risk = zeros(Num,1); 
% Annual income matrix: length (total number of houses) 
IncomeEach = zeros(Num,1); 
% Total wealth matrix: length (total number of houses) 
WealthEach = zeros(Num,1); 
% Normalized number of people in each income range: row (census tract)
% column (income range) 
NormIncome = zeros(size(Income)); 
% Expected losses of each building type given wind intensity:
% row (building type) column (wind speed)
EL = zeros(size(Total_EL)); 
% Normalized expected losses of each building type given wind intensity:
% row (building type) column (wind speed)
Norm_EL = zeros(size(Total_EL)); 
% Mean values of loss functions: length (building types)
mu_Loss = zeros(size(TypeBuilding,2),1); 
% Standard deviation of loss functions: length (building types) 
sigma_Loss = zeros(size(TypeBuilding,2),1); 
% Deductibles: length (total number of houses) 
Deduct = zeros(Num,1); 
% Index matrix showing the purchase of insurance: length (total number of
% houses)
Index_I = zeros(Num,1);
% Expected utilities 1 Equation (1a): row (total number of houses) column
% (number of intervals of coverage)
EU1 = zeros(Num,N_int_C); 
% Expected utilities 2 Equation (1b): length (total number of houses) 
EU2 = zeros(Num,1); 
% Index of insurance coverage selected: length (total number of houses)
Max_index = zeros(Num,1); 
% Inusrance coverage selected: length (total number of houses)
Selected_C = zeros(Num,1); 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 1. Hazard Curve 
% 1.1. Find the parameters of Weibull distribution (see Ref 1. Li and Ellingwood,
% 2006)
y = zeros(size(T)); 
for i = 1:length(T)
    y(i) = log(-log(1/T(i))); 
end

log_v = log(v); 
P = polyfit(log_v,y,1);
alpha = P(1); u = exp(-P(2)/alpha); 

int_v = linspace(0.1,200,200); 
y_axis_hazard = zeros(size(int_v)); 
for i = 1:length(y_axis_hazard)
    y_axis_hazard(i) = 1-wblcdf(int_v(i),u,alpha);
end
% 1.2. Find the annual mean rate of occurrence of storms in the study region
% Based on HAZUS-MH damage and loss functions, wind speed greater than 60
% m/h induces the damage and loss of a property. Thus, the probabiltiy of
% exceedance at v = 60 mph is assumed to be the mean rate of occurrence in
% this study
annual_p = y_axis_hazard(find(abs(int_v-60)<0.5)); 

% 2. Assign census tract and building type for each building 
numcount = 0; 
for i = 1:size(TypeBuilding,1)
    for j = 1:size(TypeBuilding,2)        
        for k = numcount+1:TypeBuilding(i,j)+numcount 
            IndexBldg(k,1) = i;
            IndexBldg(k,2) = j;             
        end
        numcount = TypeBuilding(i,j)+numcount; 
    end
end

% 3. Assign a threshold of concern p* to each household 
norm_p = makedist('Normal','mu',annual_p,'sigma',COV_p*annual_p); 
trunc_p = truncate(norm_p,lower_p,upper_p); % truncated normal distribution
u_p = rand(Num,1); 
for i = 1:length(p_star)
    p_star(i) = icdf(trunc_p,u_p(i)); 
end
clear u_p 

% 4. Assign alpha to each household 
norm_risk = makedist('Normal','mu',mean_risk,'sigma',cov_risk*mean_risk);  
% Truncate the distribution to have a lower limit of 0 and an upper limit of 1. 
trunc_risk = truncate(norm_risk,min_risk,max_risk);  % t2 is the distribution of alpha
u_risk = rand(Num,1);
for i=1:length(risk)
    risk(i) = icdf(trunc_risk,u_risk(i));
end
clear u_risk 

% 5. Assign individual household's income and the associated total wealth 
% Find the relationship between the annual income and total wealth 
income = [15.2 31.4 52.7 86.1 136.7]; 
wealth = [6.6 31.5 81.9 167.3 389.7];
corr_inc_wealth = polyfit(income,wealth,1); 
clear income; clear wealth; 

for i = 1:size(NormIncome,1)
    temp_cum = 0;
    for j = 1:size(NormIncome,2) 
        temp_cum = temp_cum+Income(i,j); 
        NormIncome(i,j) = temp_cum/sum(Income(i,:)); 
    end
    clear temp_cum 
end

for i = 1:length(IncomeEach)
    randn = rand(1); 
    temp_census = IndexBldg(i,1); 
    if randn <= NormIncome(temp_census,1) 
        IncomeEach(i) = IncomeRange(1)+(IncomeRange(2)-IncomeRange(1))*rand(1); 
    elseif NormIncome(temp_census,1) < randn && randn <= NormIncome(temp_census,2) 
        IncomeEach(i) = IncomeRange(2)+(IncomeRange(3)-IncomeRange(2))*rand(1);
    elseif NormIncome(temp_census,2) < randn && randn <= NormIncome(temp_census,3) 
        IncomeEach(i) = IncomeRange(3)+(IncomeRange(4)-IncomeRange(3))*rand(1);
    elseif NormIncome(temp_census,3) < randn && randn <= NormIncome(temp_census,4) 
        IncomeEach(i) = IncomeRange(4)+(IncomeRange(5)-IncomeRange(4))*rand(1);
    elseif NormIncome(temp_census,4) < randn && randn <= NormIncome(temp_census,5) 
        IncomeEach(i) = IncomeRange(5)+(IncomeRange(6)-IncomeRange(5))*rand(1);        
    elseif NormIncome(temp_census,5) < randn && randn <= NormIncome(temp_census,6) 
        IncomeEach(i) = IncomeRange(6)+(IncomeRange(7)-IncomeRange(6))*rand(1); 
    elseif NormIncome(temp_census,6) < randn && randn <= NormIncome(temp_census,7) 
        IncomeEach(i) = IncomeRange(7)+(IncomeRange(8)-IncomeRange(7))*rand(1);       
    elseif NormIncome(temp_census,7) < randn && randn <= NormIncome(temp_census,8) 
        IncomeEach(i) = IncomeRange(8)+(IncomeRange(9)-IncomeRange(8))*rand(1);  
    elseif NormIncome(temp_census,8) < randn && randn <= NormIncome(temp_census,9) 
        IncomeEach(i) = IncomeRange(9)+(IncomeRange(10)-IncomeRange(9))*rand(1);
    end
    IncomeEach(i) = IncomeEach(i)*Inf_2018/Inf_2010; 
    WealthEach(i) = corr_inc_wealth(1)*IncomeEach(i)+corr_inc_wealth(2); 
    clear randn
    clear temp_census
end

% 6. Loss functions
% Expected losses of each type of building 
for i = 1:size(EL,1)
    for j = 1:size(EL,2)
        EL(i,j) = Total_EL(i,j)/sum(TypeBuilding(:,i));
    end
end

for i = 1:size(Norm_EL,1)
    for j = 1:size(Norm_EL,2)
        Norm_EL(i,j) = Total_EL(i,j)/sum(TotalBldgValue(:,i)); 
    end
end

for i = 1:length(mu_Loss)
        P = polyfit(norminv(Norm_EL(i,:)),log(v),1);
        sigma_Loss(i) = P(1); 
        mu_Loss(i) = P(2); 
end

%{
x = 60:1:250;
y = zeros(size(Norm_EL,2),length(x));
for i = 1:size(Norm_EL,2)
    y(i,:) = logncdf(x,mu_Loss(i),sigma_Loss(i)); 
end

figure(1)
plot(x,y(1,:)); hold on; plot(v,Norm_EL(1,:),'ro'); hold off
figure(2)
plot(x,y(2,:)); hold on; plot(v,Norm_EL(2,:),'ro'); hold off
figure(3)
plot(x,y(3,:)); hold on; plot(v,Norm_EL(3,:),'ro'); hold off
figure(4)
plot(x,y(4,:)); hold on; plot(v,Norm_EL(4,:),'ro'); hold off
figure(5)
plot(x,y(5,:)); hold on; plot(v,Norm_EL(5,:),'ro'); hold off
figure(6)
plot(x,y(6,:)); hold on; plot(v,Norm_EL(6,:),'ro'); hold off
figure(7)
plot(x,y(7,:)); hold on; plot(v,Norm_EL(7,:),'ro'); hold off

figure(8)
plot(x,y(1,:),'b-'); hold on; plot(x,y(2,:),'r-'); hold on; plot(x,y(3,:),'c-'); hold on; plot(x,y(4,:),'k-'); plot(x,y(5,:),'m-'); plot(x,y(6,:),'k--'); plot(x,y(7,:),'b--');hold off;
grid on
%}

% 7. Deductible 
for i = 1:length(Deduct)
    Deduct(i) = BldgValue(IndexBldg(i,1),IndexBldg(i,2))*D; 
end

% 8. Utility model (no used)
hazard = linspace(min_H,max_H,N_int_H);
for k = 1:Num
    coverage = linspace(min_C,BldgValue(IndexBldg(k,1),IndexBldg(k,2)),N_int_C);

        for j = 1:length(coverage)
            Fterm1 = 0; % first term in Equation (1a)
            Fterm2 = 0; % first term in Equation (1b)
            for i = 1:length(hazard)
                if i == length(hazard)
                    temp_v = 200;
                    temp_p = 1-wblcdf(temp_v,u,alpha);
                else
                    temp_v = (hazard(i)+hazard(i+1))/2;
                    temp_p = wblcdf(hazard(i+1),u,alpha)-wblcdf(hazard(i),u,alpha);
                end
                temp_type = IndexBldg(k,2);
                temp_EL = logncdf(temp_v,mu_Loss(temp_type),sigma_Loss(temp_type))*BldgValue(IndexBldg(k,1),IndexBldg(k,2));
                % Claim payment 
                if temp_EL <= Deduct(k)
                    temp_CP = 0;
                elseif temp_EL > coverage(j)
                    temp_CP = coverage(j)-Deduct(k);
                else
                    temp_CP = temp_EL-Deduct(k);
                end
                temp_GA1 = 0; temp_GA2 = 0; 
                %{
                % Government assistance
                if hazard(i) <= 80
                    temp_GA1 = 0;
                else
                    if temp_EL-temp_CP < GA_max
                        temp_GA1 = temp_EL-temp_CP;
                    else
                        temp_GA1 = GA_max;
                    end
                end
                if hazard(i) <= 80
                    temp_GA2 = 0;
                else
                    if temp_EL < GA_max
                        temp_GA2 = temp_EL;
                    else
                        temp_GA2 = GA_max;
                    end
                end
                %}
                % Fterm1 = Fterm1+temp_p*(WealthEach(k)-temp_EL+temp_CP-z*coverage(j)+temp_GA1);
                % Fterm2 = Fterm2+temp_p*(WealthEach(k)-temp_EL+temp_GA2);
                % Utility function applied
                Fterm1 = Fterm1+temp_p*((WealthEach(k)-temp_EL+temp_CP-z*coverage(j)+temp_GA1).^risk(k));
                Fterm2 = Fterm2+temp_p*((WealthEach(k)-temp_EL+temp_GA2).^risk(k));
                clear temp_p
                clear temp_v
                clear temp_type
                clear temp_EL
                clear temp_CP
                clear temp_GA1
                clear temp_GA2
            end
            %Sterm1 = (1-annual_p)*(WealthEach(k)-z*coverage(j));            
            % Utility function applied
            Sterm1 = (1-annual_p)*((WealthEach(k)-z*coverage(j)).^risk(k));  
            EU1(k,j) = Fterm1+Sterm1;        
        end
        Max_index(k) = find(EU1(k,:)==max(EU1(k,:)));
        Selected_C(k) = max(EU1(k,:));
        % Sterm2 = (1-annual_p)*WealthEach(k);
        % Utility function applied
        Sterm2 = (1-annual_p)*(WealthEach(k).^risk(k));
        EU2(k) = Fterm2+Sterm2;
       
   
    clear coverage 
end

% Number of buildings in each type (WSF1, WSF2, WMUH1, WMUH2, MSF1, MSF2,
% CERBL): row (census tract) column (building type) 
save TypeBuilding
% Index matrix: row(number of bldgs) column (building type) the first column
% indicates census tract and the second column indicates building type
save IndexBldg
% Building value for each type of building in each census tract: row
% (census tract) column (building type)
save BldgValue
% Deductibles: length (total number of houses) 
save Deduct  

% Total wealth matrix: length (total number of houses) 
save WealthEach  
save Selected_C
