clear all; close all; clc 

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

% Paper: journal paper for catastrophe risk insurance 
% Case study: windstorm insurance in miami county, FL

% Phase 2: Ex-Post Community-Level Recovery Model  

% NOTE: All monetary units are thousand dollars. 
% NOTE: This code should be run after running Phase1.m,
% DamageLossRelationship.m, and damageCell.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% INPUT PARAMETERS

% Repair Time matrix: adapted from Table 3 in the paper 
% RT_ratio: building loss ratio 
% RT_time: first row (repair time for single family dwelling) second row
% (repair time for multi family dwelling) 
RT_ratio = [0 0.02 0.10 0.50 1]; 
RT_time = [0 2 30 90 180; 0 5 30 120 240]; 


% Delay Time 
% Inspection Time matrix: row(2 DS: [DS0 DS1] and [DS2 DS3 DS4]) column
% (median and cov of lognormal distributions of Inspection)
Insp_time = [0 0; 5*7 0.54]; % units: weeks
% Engineering Mobilization and Review/Redesign Time matrix: row(4 DS:
% [DS0], [DS1], [DS2 DS3], and [DS4]) column (median and cov of lognormal
% distributions of Inspection)
Eng_time = [0 0; 6*7 0.40; 12*7 0.4; 50*7 0.32]; % units: weeks
% Financing Time matrix: row (insurance and government assistance) column
% (median and cov of lognormal distributions)
Fin_time = [6*7 1.11; 48*7 0.57]; 
% Contractor Mobilization Time matrix: row (3 DS: [DS0], [DS1 DS2], and
% [DS3 DS4]) column (median and cov of lognormal distributions)
Contr_time = [0 0; 11*7 0.43; 23*7 0.41]; 
% Permitting Time matrix: row (3 DS: [DS0], [DS1 DS2], and [DS3 DS4])
% column (median and cov of lognormal distributions)
Perm_time = [0 0; 1*7 0.86; 8*7 0.32]; 


% Repair Time matrix: row(number of building type) column (median and cov
% of lognormal distributions of single family building repair time; median
% and cov of lognormal distributions of multi family building repair time)

% Damage cell: saves the percentage of building in each damage state
% following hypothetic  Hurricane 
% length (the number of building type) 
% in each cell, there is a matrix: row (census tract) column (percentage of
% buildings in each damage state)
load Cell_Damage 
% Cell_DSEL_para contains the lognormal distribution parameters of each
% damage state of each building type in each census tract 
% row (census tract) column (building type)
% each cell has a matrix: row (damage states) column (lognormal
% distribution parametsr: median and standard deviation)
load Cell_DSEL_para 
% Number of buildings in each type (WSF1, WSF2, WMUH1, WMUH2, MSF1, MSF2,
% CERBL): row (census tract) column (building type) 
load TypeBuilding
% Index matrix: row(number of bldgs) column (building type) the first column
% indicates census tract and the second column indicates building type
load IndexBldg
% Building value for each type of building in each census tract: row
% (census tract) column (building type)
load BldgValue
% Deductibles: length (total number of houses) 
load Deduct  
% Insurance coverage: length (total number of houses)
load Selected _C
% Total wealth matrix: length (total number of houses) 
load WealthEach  

min_C=6.423;

t0 = 10; % the time of occurrence of hypothetic Hurricane 
t_year = 3; % year considered 
t_days = t_year*365; % days considered 
T_inf = t_days+1; % if a homeowner files for bankruptcy, he/she moves out of 
               % town and its financing time is treated as an infinited number 
Num = sum(TypeBuilding(:)); 
GA_max = 33.792; % Maximum government assistance
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ZERO MATRICES 

% Cumulative Cell Damage matrix: cumulative percentage of building in each
% damage state following hypothetic Hurricane 
% length (the number of building type) 
% in each cell, there is a matrix: row (census tract) column (percentage of
% buildings in each damage state)
Cell_CumDam = cell(length(Cell_Damage),1); 
% Damage state matrix: save the damage state of each building following the
% scenario hypothetic hurricane
% row (total number of buildings) column (number of days considered) 
DS = zeros(Num,t_days); 
% Actual economic loss vector: compute economic losses of each building using
% the relationship between damage state and economic losses and save it (in
% thousdand dollars)
EL_act = zeros(Num,1); 
% Repair time vector: compute the repair time which is as a function of
% building loss ratio (HAZUS-MH Technical Manual) 
% length (the total number of building) 
RT_total = zeros(Num,1); 
% Inspection delay: compute the delay time of each building
T_insp = zeros(Num,1); 
% Financing delay: compute the delay time of each building 
T_fin = zeros(Num,1); 
% Engineering mobilization and review/redesign time of each building 
T_eng = zeros(Num,1);
% Permitting time of each building 
T_perm = zeros(Num,1); 
% Contractor mobilization of each building
T_contr = zeros(Num,1); 
% Bankruptcy matrix: record if the homeowner files for bankruptcy
Bankruptcy = zeros(Num,1); 
% Total delay time matrix: find the maximum time induced by three different
% delay sequences
DT_total = zeros(Num,1); 
% Recovery matrix: record if the building is fully repaired (index = 1) or
% not (index = 0); will be used to construct community-level post-disaster
% recovery curve 
% row (total number of buildings) column (number of days considered) 
Index_recovery = zeros(Num,t_days); 
% Community recovery matrix: shows the number of house buildings that have
% DM_0
comm_recovery = zeros(t_days,1); 
%claim
claim=zeros(Num,1);

% Save random number generated for the purpose of validation 
randk = zeros(Num,1);
randm = zeros(Num,1); 

% only consider the wealth greater than 0
for i=1:Num
    if WealthEach(i)<0
        WealthEach(i)=0;
    else
        WealthEach(i)= WealthEach(i);
    end
end

% 1. Determine the damage state of each building following the scenario
% hypothetic hurricane 

for i = 1:length(Cell_Damage)
    Cell_CumDam{i} = cumsum(Cell_Damage{i},2);
end
for i = 1:Num
    DS(i,1:t0) = 0;
    Index_recovery (i,1:t0) = 1;
    rand_num = rand(1);
    if rand_num <= Cell_CumDam{IndexBldg(i,2)}(IndexBldg(i,1),1)
        DS(i,t0+1) = 0; % no damage
        Index_recovery(i,t0+1) = 1;
    elseif Cell_CumDam{IndexBldg(i,2)}(IndexBldg(i,1),1) < rand_num && rand_num <= Cell_CumDam{IndexBldg(i,2)}(IndexBldg(i,1),2)
        DS(i,t0+1) = 1; % minor damage
        Index_recovery(i,t0+1) = 0;
    elseif Cell_CumDam{IndexBldg(i,2)}(IndexBldg(i,1),2) < rand_num && rand_num <= Cell_CumDam{IndexBldg(i,2)}(IndexBldg(i,1),3)
        DS(i,t0+1) = 2; % moderate damage
        Index_recovery(i,t0+1) = 0;
    elseif Cell_CumDam{IndexBldg(i,2)}(IndexBldg(i,1),3) < rand_num && rand_num <= Cell_CumDam{IndexBldg(i,2)}(IndexBldg(i,1),4)
        DS(i,t0+1) = 3; % severe damage
        Index_recovery(i,t0+1) = 0;
    else
        DS(i,t0+1) = 4; % destruction
        Index_recovery(i,t0+1) = 0;
    end
    clear rand_num
end

% 2. Compute the expected losses of each building using the damage-loss
% relationship
for i = 1:Num
    temp_tract = IndexBldg(i,1);
    temp_type = IndexBldg(i,2);
    temp_DS = DS(i,t0+1);
    if temp_DS == 0
        EL_act(i) = 0;
    elseif temp_DS == 4
        EL_act(i) = BldgValue(temp_tract,temp_type);
    else
        temp_mean = Cell_DSEL_para{temp_tract,temp_type}(temp_DS+1,1);
        temp_std = Cell_DSEL_para{temp_tract,temp_type}(temp_DS+1,2);
        randk(i) = rand(1,1); % for the purpose of verification
        EL_act(i) = logninv(randk(i),temp_mean,temp_std);
        if EL_act(i) > BldgValue(temp_tract,temp_type)
            EL_act(i) = BldgValue(temp_tract,temp_type);
        end
        clear temp_mean
        clear temp_std
    end
    clear temp_tract
    clear temp_type
    clear temp_DS
end

% 3. Compute repair time for each building
for i = 1:Num
    if EL_act(i) == 0
        RT_total(i) = 0;
    else
        temp_tract = IndexBldg(i,1);
        temp_type = IndexBldg(i,2);
        temp_loss_ratio = EL_act(i)/BldgValue(temp_tract,temp_type);
        if temp_type == 3 || temp_type == 4
            RT_total(i) = interp1(RT_ratio,RT_time(2,:),temp_loss_ratio);
        else
            RT_total(i) = interp1(RT_ratio,RT_time(1,:),temp_loss_ratio);
        end
    end
    clear temp_tract
    clear temp_type
    clear temp_loss_ratio
end

% 4. Compute delay time for each building

for i = 1:Num
    % Delay due to inspection
    if DS(i,t0+1) == 0 || DS(i,t0+1) == 1
        T_insp(i) = 0;
    else
        temp_mean = log(Insp_time(2,1));
        temp_std = sqrt(log((Insp_time(2,2)^2)+1));
        randm(i) = rand(1,1);
        T_insp(i) = logninv(randm(i),temp_mean,temp_std);
        clear temp_mean
        clear temp_std
    end
    % Delay due to engineering mobilization and review/redesign
    if DS(i,t0+1) == 0
        T_eng(i) = 0;
    elseif DS(i,t0+1) == 1
        temp_mean = log(Eng_time(2,1));
        temp_std = sqrt(log((Eng_time(2,2)^2)+1));
        rand_num = rand(1,1);
        T_eng(i) = logninv(rand_num,temp_mean,temp_std);
    elseif DS(i,t0+1) == 4
        temp_mean = log(Eng_time(4,1));
        temp_std = sqrt(log((Eng_time(4,2)^2)+1));
        rand_num = rand(1,1);
        T_eng(i) = logninv(rand_num,temp_mean,temp_std);
    else
        temp_mean = log(Eng_time(3,1));
        temp_std = sqrt(log((Eng_time(3,2)^2)+1));
        rand_num = rand(1,1);
        T_eng(i) = logninv(rand_num,temp_mean,temp_std);
    end
    clear temp_mean
    clear temp_std
    clear rand_num
    % Delay due to permitting
    if DS(i,t0+1) == 0
        T_perm(i) = 0;
    elseif DS(i,t0+1) == 1 || DS(i,t0+1) == 2
        temp_mean = log(Perm_time(2,1));
        temp_std = sqrt(log((Perm_time(2,2)^2)+1));
        rand_num = rand(1,1);
        T_perm(i) = logninv(rand_num,temp_mean,temp_std);
    else
        temp_mean = log(Perm_time(3,1));
        temp_std = sqrt(log((Perm_time(3,2)^2)+1));
        rand_num = rand(1,1);
        T_perm(i) = logninv(rand_num,temp_mean,temp_std);
    end
    clear temp_mean
    clear temp_std
    clear rand_num
    % Delay due to contractor mobilization
    if DS(i,t0+1) == 0
        T_contr(i) = 0;
    elseif DS(i,t0+1) == 1 || DS(i,t0+1) == 2
        temp_mean = log(Contr_time(2,1));
        temp_std = sqrt(log((Contr_time(2,2)^2)+1));
        rand_num = rand(1,1);
        T_contr(i) = logninv(rand_num,temp_mean,temp_std);
    else
        temp_mean = log(Contr_time(3,1));
        temp_std = sqrt(log((Contr_time(3,2)^2)+1));
        rand_num = rand(1,1);
        T_contr(i) = logninv(rand_num,temp_mean,temp_std);
    end
    clear temp_mean
    clear temp_std
    clear rand_num
    % Delay due to financing
    if DS(i,t0+1) == 0
        T_fin(i) = 0;
    else
        if EL_act(i) <= WealthEach(i)
            T_fin(i) = 0;
        else
           claim(i)=Selected_C(i);
                if EL_act(i) <= Deduct(i)
                    temp_CP = 0;
                elseif EL_act(i) >  claim(i)
                    temp_CP =  claim(i)-Deduct(i);
                else
                    temp_CP = EL_act(i)-Deduct(i);
                end
                % Government assistance
                if EL_act(i)-temp_CP < GA_max
                    temp_GA = EL_act(i)-temp_CP;
                else
                    temp_GA = GA_max;
                end
                % Check if the homeowner has sufficient financing to cover the
                % losses
                if EL_act(i)<=WealthEach(i)+temp_CP
                    temp_mean = log(Fin_time(1,1));
                    temp_std = sqrt(log((Fin_time(1,2)^2)+1));
                    rand_num = rand(1,1);
                    T_fin(i) = logninv(rand_num,temp_mean,temp_std);
                elseif EL_act(i) <= WealthEach(i)+temp_CP+temp_GA
                        temp_mean = log(Fin_time(2,1));
                        temp_std = sqrt(log((Fin_time(2,2)^2)+1));
                        rand_num = rand(1,1);
                        T_fin(i) = logninv(rand_num,temp_mean,temp_std);
                else
                    Bankruptcy(i) = 1;
                    T_fin(i) = T_inf;
                   
                end
                clear temp_CP
                clear temp_GA
                clear temp_mean
                clear temp_std
                clear rand_num
      
        end
    end
    DT_total(i) = T_insp(i)+max([T_fin(i),T_eng(i)+T_perm(i),T_contr(i)]);
end

% 5. Aggregate building recovery process to construct community recovery
% process

for i = 1:Num
    if Index_recovery(i,t0+1) ==0 
        temp_total_time = floor(DT_total(i)+RT_total(i)); 
        Index_recovery(i,t0+2:t0+1+temp_total_time) = 0; 
        Index_recovery(i,t0+1+temp_total_time+1:t_days) = 1; 
    else
        Index_recovery(i,t0+2:t_days) = 1; 
    end
    clear temp_total_time 
end

for i = 1:t_days
    comm_recovery(i) = sum(Index_recovery(:,i)); 
end

plot(1:1:(t_days-1),comm_recovery(1:t_days-1));
save comm_recovery


