Repository: IN-CORE/pyincore (language: python)

Issue:
Housing recovery model analysis
The housing recovery model using Census appraisal data, Galveston testbed notebook example. 
Notes: https://opensource.ncsa.illinois.edu/confluence/display/INCORE1/Housing+recovery+analysis


Discussion / hints:

Fix patch:
diff --git a/docs/source/modules.rst b/docs/source/modules.rst
index 9eed86f77..d8e82cb92 100644
--- a/docs/source/modules.rst
+++ b/docs/source/modules.rst
@@ -63,6 +63,13 @@ analyses/example
 ..  autoclass:: example.exampleanalysis.ExampleAnalysis
     :members:
 
+analyses/housingrecovery
+========================
+..  autoclass:: housingrecovery.housingrecovery.HousingRecovery
+    :members:
+..  autoclass:: housingrecovery.housingrecoveryutil.HousingRecoveryUtil
+    :members:
+
 analyses/housingrecoverysequential
 ==================================
 ..  autoclass:: housingrecoverysequential.housingrecoverysequential.HousingRecoverySequential
diff --git a/pyincore/analyses/housingrecovery/__init__.py b/pyincore/analyses/housingrecovery/__init__.py
new file mode 100644
index 000000000..89d28bc27
--- /dev/null
+++ b/pyincore/analyses/housingrecovery/__init__.py
@@ -0,0 +1,8 @@
+# Copyright (c) 2019 University of Illinois and others. All rights reserved.
+#
+# This program and the accompanying materials are made available under the
+# terms of the Mozilla Public License v2.0 which accompanies this distribution,
+# and is available at https://www.mozilla.org/en-US/MPL/2.0/
+
+
+from pyincore.analyses.housingrecovery.housingrecovery import HousingRecovery
diff --git a/pyincore/analyses/housingrecovery/housingrecovery.py b/pyincore/analyses/housingrecovery/housingrecovery.py
new file mode 100755
index 000000000..ac1687d82
--- /dev/null
+++ b/pyincore/analyses/housingrecovery/housingrecovery.py
@@ -0,0 +1,386 @@
+# This program and the accompanying materials are made available under the
+# terms of the Mozilla Public License v2.0 which accompanies this distribution,
+# and is available at https://www.mozilla.org/en-US/MPL/2.0/
+
+import sys
+import numpy as np
+import pandas as pd
+
+from pyincore import BaseAnalysis
+from pyincore.analyses.housingrecovery.housingrecoveryutil import HousingRecoveryUtil
+
+
+class HousingRecovery(BaseAnalysis):
+    """ The analysis predicts building values and value changes over time following a disaster event.
+    The model is calibrated with respect to demographics, parcel data, and building value trajectories
+    following Hurricane Ike (2008) in Galveston, Texas. The model predicts building value at the parcel
+    level for 8 years of observation. The models rely on Census (Decennial or American Community Survey, ACS)
+    and parcel data immediately prior to the disaster event (year -1) as inputs for prediction.
+
+    The Galveston, TX example makes use of 2010 Decennial Census and Galveston County Appraisal District (GCAD)
+    tax assessor data and outputs from other analysis (i.e., Building Damage, Housing Unit Allocation,
+    Population Dislocation) .
+
+    The CSV outputs of the building values for the 6 years following the disaster event (with year 0 being
+    the impact year).
+
+    Contributors
+        | Science: Wayne Day, Sarah Hamideh
+        | Implementation: Michal Ondrejcek, Santiago Núñez-Corrales, and NCSA IN-CORE Dev Team
+
+    Related publications
+        Hamideh, S., Peacock, W. G., & Van Zandt, S. (2018). Housing recovery after disasters: Primary versus
+        seasonal/vacation housing markets in coastal communities. Natural Hazards Review.
+
+    Args:
+        incore_client (IncoreClient): Service authentication.
+
+    """
+    def __init__(self, incore_client):
+        super(HousingRecovery, self).__init__(incore_client)
+
+    def get_spec(self):
+        return {
+            "name": "housing-recovery",
+            "description": "Housing Recovery Analysis",
+            "input_parameters": [
+                {
+                    "id": "base_year",
+                    "required": False,
+                    "description": "Base year is used to calculate improvement age. It needs to be set to the tax "
+                                   "assessment year representing pre-disaster building values. For example for GCAD "
+                                   "data which represents improvement valuation before Hurricane Ike impacts."
+                                   "Deafult 2008",
+                    "type": int
+                },
+                {
+                   "id": "result_name",
+                   "required": True,
+                   "description": "Result CSV dataset name",
+                   "type": str
+                }
+            ],
+            "input_datasets": [
+                {
+                   "id": "population_dislocation",
+                   "required": True,
+                   "description": "Population Dislocation aggregated to the block group level",
+                   "type": ["incore:popDislocation"]
+                },
+                {
+                   "id": "building_area",
+                   "required": True,
+                   "description": "Building square footage and damage. Damage is the actual building value loss "
+                                  "in percentage terms observed through the County Appraisal District (GCAD) data",
+                   "type": ["incore:buildingInventoryArea"]
+                },
+                {
+                   "id": "census_block_groups_data",
+                   "required": True,
+                   "description": "Census ACS data, 2010 5yr data for block groups available at IPUMS NHGIS "
+                                  "website.",
+                   "type": ["incore:censusBlockGroupsData"]
+                },
+                {
+                   "id": "census_appraisal_data",
+                   "required": True,
+                   "description": "Census data, 2010 Decennial Census District (GCAD) Census data",
+                   "type": ["incore:censusAppraisalData"]
+                }
+            ],
+            "output_datasets": [
+                {
+                   "id": "result",
+                   "description": "A csv file with the building values for the 6 years following the disaster"
+                                  "event (year -1 denotes pre-impact conditions and 0 being the impact year). "
+                                  "Index year values represent building values against a base, pre-impact value.",
+                   "type": "incore:buildingValues"
+                }
+            ]
+        }
+
+    def run(self):
+        """ Executes the Housing recovery analysis.
+
+        Returns:
+            bool: True if successful, False otherwise.
+
+        """
+        hru = HousingRecoveryUtil()
+
+        # Get base year
+        self.base_year = self.get_parameter("base_year")
+        if self.base_year is None:
+            self.base_year = hru.BASEYEAR
+        # Get desired result name
+        result_name = self.get_parameter("result_name")
+
+        # Datasets
+        pop_disl = self.get_input_dataset("population_dislocation").get_dataframe_from_csv(low_memory=False)
+        addl_structure_info = self.get_input_dataset("building_area").get_dataframe_from_csv(low_memory=False)
+        bg_mhhinc = self.get_input_dataset("census_block_groups_data").get_dataframe_from_csv(low_memory=False)
+
+        # Census data
+        vac_status = self.get_input_dataset("census_appraisal_data").get_json_reader()
+
+        # Calculate the percent vacation or seasonal housing of all housing units within a census tract
+        vac_status = self.get_vac_season_housing(vac_status)
+
+        # Get ownership by filtering vacancy codes
+        pop_disl["d_ownerocc"] = self.get_owneship(pop_disl)
+
+        # Keep only parcels with 1 building
+        building_level = pop_disl.loc[(pop_disl.bldgobs == 1)]
+        single_family = building_level.loc[(building_level.d_sf == 1)]
+
+        # Read in and clean block group level median household income
+        bg_mhhinc["mhhinck"] = bg_mhhinc["mhhinc"] / 1000
+
+        # Read in and clean additional building information, NOTE add to building inventory
+        # Create structure id for merging data
+        addl_structure_info["strctid"] = addl_structure_info["xref"].apply(lambda x: "XREF"+x)
+        hse_recov = self.merge_add_inv(single_family, addl_structure_info)
+
+        # Merge with seasonal/vacation housing Census ACS data
+        hse_recov = self.merge_seasonal_data(hse_recov, vac_status)
+
+        # Merge with BG median HH income
+        hse_recov = self.merge_block_data(hse_recov, bg_mhhinc)
+
+        # Generate minority variable
+        hse_recov["pminoritybg"] = hse_recov["phispbg"] + hse_recov["pblackbg"]
+
+        # Estimate value_loss for each parcel based on parameters from Bai, Hueste, & Gardoni (2009)
+        hse_recov = self.value_loss(hse_recov)
+
+        # Primary Housing Market (PHM)
+        hse_rec_phm = self.assemble_phm_coefs(hru, hse_recov)
+        # Seasonal/Vacation housing market (SVHM)
+        hse_rec_svhm = self.assemble_svhm_coefs(hru, hse_recov)
+
+        d_vac_np = np.tile(hse_recov["d_vacationct"], (8, 1)).T
+        # Vacation condition for all years
+
+        hse_rec_fin = np.empty(hse_recov.shape)
+        hse_rec_fin[:] = np.NaN
+        hse_rec_fin = np.where(d_vac_np == 0, np.exp(hse_rec_phm), np.NaN)
+        hse_rec_fin = np.where(d_vac_np == 1, np.exp(hse_rec_svhm), hse_rec_fin)
+
+        # Index building values against a pre-impact (-1) value
+        hse_rec_fin_index = hse_rec_fin / hse_rec_fin[:, 0:1]
+
+        # Name columns, use list of damage years (range from -1 to n).
+        bval_yr = ["bval_year_{0}".format(i) for i in hru.DMG_YEARS]
+        index_yr = ["index_{0}".format(i) for i in hru.DMG_YEARS]
+
+        hse_rec_fin1 = pd.DataFrame(hse_rec_fin, columns=bval_yr)
+        hse_rec_fin2 = pd.DataFrame(hse_rec_fin_index, columns=index_yr)
+        hse_recov = pd.concat([hse_recov, hse_rec_fin1, hse_rec_fin2], axis=1)
+
+        column_to_save = ["guid", "d_vacationct", "mhhinck", "pminoritybg", "dmg", "value_loss"] + bval_yr + index_yr
+        hse_recov[column_to_save].to_csv(result_name + ".csv", index=False, float_format="%.2f")
+
+        return True
+
+    def get_owneship(self, popd):
+        """ Filter ownership based on the vacancy codes
+        Assumption:
+        Where ownershp is "missing", let vacancy codes 0/3/4 be considered owner-occupied,
+        and 1/2/5/6/7 be considered renter-occupied.
+        It is uncertain whether vacancy codes 3,4,5,6,7 will become owner- or renter-occupied or primarily
+        one or the other.
+    .
+        Args:
+            popd (pd.DataFrame): Population dislocation results with ownership information.
+
+        Returns:
+            pd.DataFrame: Ownership data.
+
+        """
+        # Create ownership dummy variable from popd.ownership
+        own = np.where(popd["ownershp"].isnull() & popd["vacancy"] == 7, 0, np.nan)
+        own = np.where(popd["ownershp"].isnull() & popd["vacancy"] == 6, 0, own)
+        own = np.where(popd["ownershp"].isnull() & popd["vacancy"] == 5, 0, own)
+        own = np.where(popd["ownershp"].isnull() & popd["vacancy"] == 2, 0, own)
+        own = np.where(popd["ownershp"].isnull() & popd["vacancy"] == 1, 0, own)
+        own = np.where(popd["ownershp"].isnull() & popd["vacancy"] == 4, 1, own)
+        own = np.where(popd["ownershp"].isnull() & popd["vacancy"] == 3, 1, own)
+        own = np.where(popd["ownershp"].isnull() & popd["vacancy"] == 0, 1, own)
+        own = np.where(popd["ownershp"] == 2, 0, own)
+        own = np.where(popd["ownershp"] == 1, 1, own)
+
+        return own
+
+    def get_vac_season_housing(self, vacation_status):
+        """ Calculate the percent vacation or seasonal housing of all housing units within a census tract and
+        add dummy variable for census tract as a seasonal/vacation housing submarket.
+.
+        Args:
+            vacation_status (obj): Seasonal/vacation housing Census ACS data from json reader.
+
+        Returns:
+            pd.DataFrame: Seasonal/vacation housing data.
+
+        """
+        vac_status = pd.DataFrame(data=vacation_status[1:], columns=vacation_status[0])
+
+        vac_status["B25004_006E"] = vac_status["B25004_006E"].astype(int)
+        vac_status["B25004_006M"] = vac_status["B25004_006M"].astype(int)
+        vac_status["B25002_001E"] = vac_status["B25002_001E"].astype(int)
+        vac_status["B25002_001M"] = vac_status["B25002_001M"].astype(int)
+
+        # Calculate the percent vacation or seasonal housing of all housing units within a census tract
+        vac_status["pvacationct_moe"] = vac_status["B25004_006E"] / vac_status["B25002_001E"]
+        vac_status["pvacationct"] = 100 * vac_status["pvacationct_moe"]
+        vac_status["pvacationct_moe"] = vac_status["pvacationct_moe"] ** 2 * vac_status["B25002_001M"] ** 2
+        vac_status["pvacationct_moe"] = vac_status["B25004_006M"] ** 2 - vac_status["pvacationct_moe"]
+        vac_status["pvacationct_moe"] = 100 * (1 / vac_status["B25002_001E"]) * vac_status["pvacationct_moe"] ** 0.5
+
+        # dummy variable for census tract as a seasonal/vacation housing submarket
+        vac_status["d_vacationct"] = np.where(vac_status["pvacationct"] >= 50, 1, 0)
+
+        vac_status.sort_values(by="pvacationct", inplace=True, ascending=False)
+
+        return vac_status
+
+    def merge_add_inv(self, hse_rec, addl_struct):
+        """Merge study area and additional structure information.
+.
+        Args:
+            hse_rec (pd.DataFrame):  Area inventory.
+            addl_struct (pd.DataFrame):  Additional infrastructure inventory.
+
+        Returns:
+            pd.DataFrame: Final merge of two inventories.
+
+        """
+        hse_rec_merged = pd.merge(hse_rec, addl_struct, on="strctid", how="inner")
+        return hse_rec_merged
+
+    def merge_seasonal_data(self, hse_rec, vac_status):
+        """ Merge study area and with seasonal/vacation housing Census ACS data.
+.
+        Args:
+            hse_rec (pd.DataFrame): Area inventory.
+            vac_status (pd.DataFrame): Seasonal/vacation housing Census ACS data.
+
+        Returns:
+            pd.DataFrame: Final merge of two inventories.
+
+        """
+        hse_rec["tractid"] = hse_rec["tractid"].astype(str)
+        # Add county and state to trac to match hse_rec tracid (Galveston - 723900 to 48167723900)
+        vac_status["tractid"] = vac_status["state"].astype(str) + vac_status["county"].astype(str) + \
+                                vac_status["tract"].astype(str)
+
+        hse_rec_merged = pd.merge(hse_rec, vac_status, left_on="tractid", right_on="tractid", how='inner')
+        return hse_rec_merged
+
+    def merge_block_data(self, hse_rec, bg_mhhinc):
+        """ Merge block group level median household income.
+.
+        Args:
+            hse_rec (pd.DataFrame):  Area inventory.
+            bg_mhhinc (pd.DataFrame):  Block data.
+
+        Returns:
+            pd.DataFrame: Final merge of two inventories.
+
+        """
+        hse_rec_merged = pd.merge(hse_rec, bg_mhhinc, left_on="bgidstr", right_on="bgidstr", how="inner")
+        return hse_rec_merged
+
+    def value_loss(self, hse_rec):
+        """ Estimate value_loss for each parcel based on parameters from Bai, Hueste, & Gardoni (2009).
+    .
+        Args:
+            hse_rec (pd.DataFrame):  Area inventory.
+
+        Returns:
+            pd.DataFrame: Inventory with value losses.
+
+        """
+        hse_rec["value_loss"] = 100 * (hse_rec["DS_0"] * hse_rec["rploss_0"] +
+                                       hse_rec["DS_1"] * hse_rec["rploss_1"] +
+                                       hse_rec["DS_2"] * hse_rec["rploss_2"] +
+                                       hse_rec["DS_3"] * hse_rec["rploss_3"])
+        return hse_rec
+
+    def assemble_phm_coefs(self, hru, hse_rec):
+        """ Assemble Primary Housing Market (PHM) data for full inventory and all damage-related years.
+    .
+        Args:
+            hru (obj): Housing recovery utility.
+            hse_rec (pd.DataFrame):  Area inventory including losses.
+
+        Returns:
+            np.array: Final coefficients for all damage years.
+
+        """
+        dmg_years = np.array(hru.DMG_YEARS)
+        dmg_years_size = len(hru.DMG_YEARS)
+
+        coef_fin = np.empty((hse_rec.shape[0], dmg_years_size))
+        coef_fin[:] = hru.B_PHM_intercept + np.fromiter(hru.B_PHM_year.values(), dtype=float)
+
+        # Adjust build year year with damage years
+        yrbl_all = np.empty((hse_rec.shape[0], dmg_years_size))
+        yrbl_all[:] = self.base_year + dmg_years + 1
+        yrbuilt = hru.B_PHM_age * (yrbl_all - hse_rec["effyrbuilt"].to_numpy()[:, np.newaxis])
+
+        # Square meters, use vector (1x8) with B_PHM_sqm
+        sqmeter = np.full((1, dmg_years_size), hru.B_PHM_sqm) * hse_rec["sqmeter"].to_numpy()[:, np.newaxis]
+
+        if "dmg" in hse_rec.columns:
+            dmg_loss = np.fromiter(hru.B_PHM_dmg_year.values(), dtype=float) * \
+                       hse_rec["dmg"].to_numpy()[:, np.newaxis]
+        else:
+            dmg_loss = np.fromiter(hru.B_PHM_dmg_year.values(), dtype=float) * \
+                       hse_rec["value_loss"].to_numpy()[:, np.newaxis]
+        d_owner = np.fromiter(hru.B_PHM_own_year.values(), dtype=float) * \
+                  hse_rec["d_ownerocc"].to_numpy()[:, np.newaxis]
+        mhhinck = np.fromiter(hru.B_PHM_inc_year.values(), dtype=float) * \
+                  hse_rec["mhhinck"].to_numpy()[:, np.newaxis]
+        pminrbg = np.fromiter(hru.B_PHM_min_year.values(), dtype=float) * \
+                  hse_rec["pminoritybg"].to_numpy()[:, np.newaxis]
+
+        return coef_fin + yrbuilt + sqmeter + d_owner + dmg_loss + mhhinck + pminrbg
+
+    def assemble_svhm_coefs(self, hru, hse_rec):
+        """ Assemble Seasonal/Vacation housing market (SVHM) data for full inventory and all damage-related years.
+    .
+        Args:
+            hru (obj): Housing recovery utility.
+            hse_rec (pd.DataFrame):  Area inventory including losses.
+
+        Returns:
+            np.array: Final coefficients for all damage years.
+
+        """
+        dmg_years = np.array(hru.DMG_YEARS)
+        dmg_years_size = len(hru.DMG_YEARS)
+
+        coef_fin = np.empty((hse_rec.shape[0], dmg_years_size))
+        coef_fin[:] = hru.B_SVHM_intercept + np.fromiter(hru.B_SVHM_year.values(), dtype=float)
+
+        # Adjust build year year with damage years
+        yrbl_all = np.empty((hse_rec.shape[0], dmg_years_size))
+        yrbl_all[:] = self.base_year + dmg_years + 1
+        yrbuilt = hru.B_SVHM_age * (yrbl_all - hse_rec["effyrbuilt"].to_numpy()[:, np.newaxis])
+
+        # Square meters, use vector (1x8) with B_PHM_sqm
+        sqmeter = np.full((1, dmg_years_size), hru.B_SVHM_sqm) * hse_rec["sqmeter"].to_numpy()[:, np.newaxis]
+
+        if "dmg" in hse_rec.columns:
+            dmg_loss = np.fromiter(hru.B_SVHM_dmg_year.values(), dtype=float) * \
+                       hse_rec["dmg"].to_numpy()[:, np.newaxis]
+        else:
+            dmg_loss = np.fromiter(hru.B_SVHM_dmg_year.values(), dtype=float) * \
+                       hse_rec["value_loss"].to_numpy()[:, np.newaxis]
+        d_owner = np.fromiter(hru.B_SVHM_own_year.values(), dtype=float) * \
+                  hse_rec["d_ownerocc"].to_numpy()[:, np.newaxis]
+        mhhinck = np.fromiter(hru.B_SVHM_inc_year.values(), dtype=float) * \
+                  hse_rec["mhhinck"].to_numpy()[:, np.newaxis]
+
+        return coef_fin + yrbuilt + sqmeter + dmg_loss + d_owner + mhhinck
+
diff --git a/pyincore/analyses/housingrecovery/housingrecoveryutil.py b/pyincore/analyses/housingrecovery/housingrecoveryutil.py
new file mode 100644
index 000000000..93c1fea64
--- /dev/null
+++ b/pyincore/analyses/housingrecovery/housingrecoveryutil.py
@@ -0,0 +1,128 @@
+# This program and the accompanying materials are made available under the
+# terms of the Mozilla Public License v2.0 which accompanies this distribution,
+# and is available at https://www.mozilla.org/en-US/MPL/2.0/
+
+class HousingRecoveryUtil:
+    BASEYEAR = 2008
+
+    DMG_YEARS = [-1, 0, 1, 2, 3, 4, 5, 6]
+
+    # The models presented in this section are based on models developed in Hamideh et al. (2018).
+    # The modifications from the original Galveston models include (1) combining non-Hispanic Black
+    # and Hispanic block group percentages into a minority block group percentage variable,
+    # and (2) having both parcel and block group level intercepts. The original model only had parcel level
+    # intercepts. Considering that multiple levels of data are used (parcel and block group), it is appropriate
+    # to have both levels of intercepts.
+
+    # Primary housing market (PHM) model coefficients
+    B_PHM_intercept = 11.407594
+
+    # Year indicator dummy variables
+    B_PHM_year = {}
+    B_PHM_year[-1] = 0.000000  # year -1, tax assessment immediately before disaster
+    B_PHM_year[0] = 0.263500  # year  0, tax assessment immediately after disaster, damage year
+    B_PHM_year[1] = 0.147208  # year +1
+    B_PHM_year[2] = 0.110004  # year +2
+    B_PHM_year[3] = 0.122228  # year +3
+    B_PHM_year[4] = 0.102886  # year +4
+    B_PHM_year[5] = 0.102806  # year +5
+    B_PHM_year[6] = 0.276548  # year +6
+
+    # House Characteristics
+    B_PHM_age = -0.030273
+    B_PHM_sqm = 0.005596
+
+    # Damage and year dummy interactions
+    B_PHM_dmg_year = {}
+    B_PHM_dmg_year[-1] = 0.000000
+    B_PHM_dmg_year[0] = -0.037458  # base effect
+    B_PHM_dmg_year[1] = 0.009310 + B_PHM_dmg_year[0]
+    B_PHM_dmg_year[2] = 0.009074 + B_PHM_dmg_year[0]
+    B_PHM_dmg_year[3] = 0.010340 + B_PHM_dmg_year[0]
+    B_PHM_dmg_year[4] = 0.011293 + B_PHM_dmg_year[0]
+    B_PHM_dmg_year[5] = 0.012178 + B_PHM_dmg_year[0]
+    B_PHM_dmg_year[6] = 0.012188 + B_PHM_dmg_year[0]
+
+    # Owner-occupied and year dummy interactions
+    B_PHM_own_year = {}
+    B_PHM_own_year[-1] = 0.017153   # base effect
+    B_PHM_own_year[0] = 0.129077 + B_PHM_own_year[-1]
+    B_PHM_own_year[1] = 0.188217 + B_PHM_own_year[-1]
+    B_PHM_own_year[2] = 0.235435 + B_PHM_own_year[-1]
+    B_PHM_own_year[3] = 0.246437 + B_PHM_own_year[-1]
+    B_PHM_own_year[4] = 0.261220 + B_PHM_own_year[-1]
+    B_PHM_own_year[5] = 0.281015 + B_PHM_own_year[-1]
+    B_PHM_own_year[6] = 0.265465 + B_PHM_own_year[-1]
+
+    # Median household income and year dummy interactions
+    B_PHM_inc_year = {}
+    B_PHM_inc_year[-1] = 0.002724   # base effect
+    B_PHM_inc_year[0] = 0.001190 + B_PHM_inc_year[-1]
+    B_PHM_inc_year[1] = 0.001480 + B_PHM_inc_year[-1]
+    B_PHM_inc_year[2] = 0.001746 + B_PHM_inc_year[-1]
+    B_PHM_inc_year[3] = 0.001494 + B_PHM_inc_year[-1]
+    B_PHM_inc_year[4] = 0.001480 + B_PHM_inc_year[-1]
+    B_PHM_inc_year[5] = 0.001849 + B_PHM_inc_year[-1]
+    B_PHM_inc_year[6] = 0.001400 + B_PHM_inc_year[-1]
+
+    # Block Group percent Minority and year dummy interactions
+    B_PHM_min_year = {}
+    B_PHM_min_year[-1] = -0.004783   # base effect
+    B_PHM_min_year[0] = 0.005609 + B_PHM_min_year[-1]
+    B_PHM_min_year[1] = 0.007343 + B_PHM_min_year[-1]
+    B_PHM_min_year[2] = 0.007459 + B_PHM_min_year[-1]
+    B_PHM_min_year[3] = 0.006421 + B_PHM_min_year[-1]
+    B_PHM_min_year[4] = 0.006308 + B_PHM_min_year[-1]
+    B_PHM_min_year[5] = 0.005770 + B_PHM_min_year[-1]
+    B_PHM_min_year[6] = 0.005620 + B_PHM_min_year[-1]
+
+    # Seasonal/Vacation Housing Market (SVHM) model coefficients
+    B_SVHM_intercept = 11.125636
+
+    # Year indicator dummy variables
+    B_SVHM_year =  {}
+    B_SVHM_year[-1] = 0.000000  # year -1, tax assessment immediately before disaster
+    B_SVHM_year[0] = 1.489008  # year  0, tax assessment immediately after disaster, damage year
+    B_SVHM_year[1] = 1.858770  # year +1
+    B_SVHM_year[2] = 2.163492  # year +2
+    B_SVHM_year[3] = 2.071690  # year +3
+    B_SVHM_year[4] = 2.110245  # year +4
+    B_SVHM_year[5] = 2.208577  # year +5
+    B_SVHM_year[6] = 2.273867  # year +6
+
+    # House Characteristics
+    B_SVHM_age = -0.014260
+    B_SVHM_sqm = 0.005505
+
+    # Damage and year dummy interactions
+    B_SVHM_dmg_year = {}
+    B_SVHM_dmg_year[-1] = 0.000000
+    B_SVHM_dmg_year[0] = -0.051913  # base effect
+    B_SVHM_dmg_year[1] = 0.001504 + B_SVHM_dmg_year[0]
+    B_SVHM_dmg_year[2] = -0.004172 + B_SVHM_dmg_year[0]
+    B_SVHM_dmg_year[3] = -0.002016 + B_SVHM_dmg_year[0]
+    B_SVHM_dmg_year[4] = 0.001687 + B_SVHM_dmg_year[0]
+    B_SVHM_dmg_year[5] = 0.000876 + B_SVHM_dmg_year[0]
+    B_SVHM_dmg_year[6] = 0.001129 + B_SVHM_dmg_year[0]
+
+    # Owner-occupied and year dummy interactions
+    B_SVHM_own_year = {}
+    B_SVHM_own_year[-1] = -0.017167   # base effect
+    B_SVHM_own_year[0] = 0.043263 + B_SVHM_own_year[-1]
+    B_SVHM_own_year[1] = 0.003315 + B_SVHM_own_year[-1]
+    B_SVHM_own_year[2] = 0.034372 + B_SVHM_own_year[-1]
+    B_SVHM_own_year[3] = -0.014929 + B_SVHM_own_year[-1]
+    B_SVHM_own_year[4] = -0.021078 + B_SVHM_own_year[-1]
+    B_SVHM_own_year[5] = 0.001161 + B_SVHM_own_year[-1]
+    B_SVHM_own_year[6] = 0.017562 + B_SVHM_own_year[-1]
+
+    # Median household income and year dummy interactions
+    B_SVHM_inc_year = {}
+    B_SVHM_inc_year[-1] = 0.003786   # base effect
+    B_SVHM_inc_year[0] = -0.013662 + B_SVHM_inc_year[-1]
+    B_SVHM_inc_year[1] = -0.017401 + B_SVHM_inc_year[-1]
+    B_SVHM_inc_year[2] = -0.021541 + B_SVHM_inc_year[-1]
+    B_SVHM_inc_year[3] = -0.019756 + B_SVHM_inc_year[-1]
+    B_SVHM_inc_year[4] = -0.020789 + B_SVHM_inc_year[-1]
+    B_SVHM_inc_year[5] = -0.021555 + B_SVHM_inc_year[-1]
+    B_SVHM_inc_year[6] = -0.019781 + B_SVHM_inc_year[-1]


Test patch:
diff --git a/tests/pyincore/analyses/housingrecovery/test_housingrecovery.py b/tests/pyincore/analyses/housingrecovery/test_housingrecovery.py
new file mode 100644
index 000000000..7eaf767e0
--- /dev/null
+++ b/tests/pyincore/analyses/housingrecovery/test_housingrecovery.py
@@ -0,0 +1,182 @@
+from pyincore import IncoreClient
+from pyincore import MappingSet, FragilityService
+from pyincore.analyses.buildingdamage.buildingdamage import BuildingDamage
+from pyincore.analyses.housingunitallocation import HousingUnitAllocation
+from pyincore.analyses.populationdislocation import PopulationDislocation
+from pyincore.analyses.housingrecovery.housingrecovery import HousingRecovery
+import pyincore.globals as pyglobals
+
+
+def run_with_base_class(chained):
+    client = IncoreClient(pyglobals.INCORE_API_DEV_URL)
+
+    if chained:
+        # *************************************************************
+        # Set up and run Building damage, Galveston, TX analysis
+        # *************************************************************
+        hazard_type = "hurricane"
+        # Galveston deterministic Hurricane, 3 datasets - Kriging
+        # hazard_id = "5fa5a228b6429615aeea4410"  # prod
+        hazard_id = "5fa472033c1f0c73fe81461a"
+
+        # Galveston, TX Building inventory
+        # bldg_inv_id = "60354b6c123b4036e6837ef7"  # prod 19k buildings with age_group
+        bldg_inv_id = "602eba8bb1db9c28aef01358"
+
+        # Hurricane building mapping (with equation)
+        mapping_id = "602c381a1d85547cdc9f0675" # dev & prod
+        fragility_service = FragilityService(client)
+        mapping_set = MappingSet(fragility_service.get_mapping(mapping_id))
+
+        bldg_dmg = BuildingDamage(client)
+        bldg_dmg.load_remote_input_dataset("buildings", bldg_inv_id)
+        bldg_dmg.set_input_dataset("dfr3_mapping_set", mapping_set)
+
+        result_name = "Galveston_bldg_dmg_result"
+
+        bldg_dmg.set_parameter("fragility_key", "Hurricane SurgeLevel and WaveHeight Fragility ID Code")
+        bldg_dmg.set_parameter("result_name", result_name)
+        bldg_dmg.set_parameter("hazard_type", hazard_type)
+        bldg_dmg.set_parameter("hazard_id", hazard_id)
+        bldg_dmg.set_parameter("num_cpu", 4)
+
+        # Run Building damage
+        bldg_dmg.run_analysis()
+
+        # Retrieve result dataset
+        building_dmg_result = bldg_dmg.get_output_dataset("ds_result")
+
+        # Convert dataset to Pandas DataFrame
+        # bdmg_df = building_dmg_result.get_dataframe_from_csv(low_memory=False)
+        # print(bdmg_df[["guid", "DS_0", "DS_1", "DS_2", "DS_3", "haz_expose"]])
+
+        # *************************************************************
+        # Set up and run Housing Unit Allocation (HUA) analysis
+        # *************************************************************
+        # Housing unit inventory
+        # housing_unit_inv = "5fc6ab1cd2066956f49e7a03" # prod
+        housing_unit_inv = "5fc6a757eba0eb743dc7739f"
+        # Address point inventory
+        # address_point_inv = "5fc6aadcc38a0722f563392e"
+        address_point_inv = "5fc6a4ddeba0eb743dc77301"
+
+        # Create housing allocation
+        hua = HousingUnitAllocation(client)
+
+        # Load input dataset
+        hua.load_remote_input_dataset("housing_unit_inventory", housing_unit_inv)
+        hua.load_remote_input_dataset("address_point_inventory", address_point_inv)
+        hua.load_remote_input_dataset("buildings", bldg_inv_id)
+
+        # Specify the result name
+        result_name = "Galveston_HUA"
+
+        seed = 1238
+        iterations = 1
+
+        # Set analysis parameters
+        hua.set_parameter("result_name", result_name)
+        hua.set_parameter("seed", seed)
+        hua.set_parameter("iterations", iterations)
+
+        # Run Housing unit allocation analysis
+        hua.run_analysis()
+
+        # Retrieve result dataset
+        hua_result = hua.get_output_dataset("result")
+
+        # Convert dataset to Pandas DataFrame
+        # hua_df = hua_result.get_dataframe_from_csv(low_memory=False)
+        # Display top 5 rows of output data
+        # print(hua_df[["guid", "numprec", "incomegroup"]].head())
+
+        # *************************************************************
+        # Set up and run Population Dislocation analysis
+        # *************************************************************
+
+        # Block group data, IN-CORE_BGMAP_2021-01-19_GalvestonTX
+        # bg_data = "603545f2dcda03378087e708" # prod
+        bg_data = "6035445b1e456929c86094c8"
+        # Value loss parameters DS 0-3
+        # value_loss = "60354810e379f22e16560dbd"
+        value_loss = "613901d9ca423a005e422feb"
+
+        pop_dis = PopulationDislocation(client)
+
+        pop_dis.load_remote_input_dataset("block_group_data", bg_data)
+        pop_dis.load_remote_input_dataset("value_loss_param", value_loss)
+
+        pop_dis.set_input_dataset("building_dmg", building_dmg_result)
+        pop_dis.set_input_dataset("housing_unit_allocation", hua_result)
+
+        result_name = "galveston-pop-disl-results"
+        seed = 1111
+
+        pop_dis.set_parameter("result_name", result_name)
+        pop_dis.set_parameter("seed", seed)
+
+        pop_dis.run_analysis()
+
+        # Retrieve result dataset
+        pd_result = pop_dis.get_output_dataset("result")
+
+        # Convert dataset to Pandas DataFrame
+        # pd_df = pd_result.get_dataframe_from_csv(low_memory=False)
+        # Display top 5 rows of output data
+        # print(pd_df[["guid", "dislocated"]].head())
+
+        # *************************************************************
+        # Set up and run Housing Recovery analysis
+        # *************************************************************
+
+        # Additional inventory data (assesed damage, square area)
+        bldg_sqft_id = "62193c19a42a3e546ae2d9f8"
+        # Census block groups data (required)
+        census_bg_id = "62193b7ca42a3e546ae2d9f2"
+
+        # Census appraisal file; id of external Census json is required
+        census_appr_id = "6241fbd153302c512d685181" # dev
+        result_name = "Galveston_building_values_chained"
+
+        housing_rec = HousingRecovery(client)
+
+        housing_rec.set_input_dataset("population_dislocation", pd_result)
+        housing_rec.load_remote_input_dataset("building_area", bldg_sqft_id)
+        housing_rec.load_remote_input_dataset("census_block_groups_data", census_bg_id)
+        # Census appraisal data
+        housing_rec.load_remote_input_dataset("census_appraisal_data", census_appr_id)
+
+        housing_rec.set_parameter("result_name", result_name)
+
+        housing_rec.run_analysis()
+
+    else:
+        # Run Housing Recovery analysis without chaining
+
+        # Additional inventory data (assesed damage, square area)
+        bldg_add_id = "62193c19a42a3e546ae2d9f8"
+        # Census block groups data (required)
+        census_bg_id = "62193b7ca42a3e546ae2d9f2"
+
+        # Census appraisal file; id of external Census json is required if not fips for API request
+        census_appr_id = "6241fbd153302c512d685181"
+        result_name = "Galveston_building_values"
+
+        housing_rec = HousingRecovery(client)
+
+        pop_disl_id = "623d1e1ca42a3e546aeba25f"  # dev
+        housing_rec.load_remote_input_dataset("population_dislocation", pop_disl_id)
+        housing_rec.load_remote_input_dataset("building_area", bldg_add_id)
+        housing_rec.load_remote_input_dataset("census_block_groups_data", census_bg_id)
+        # Census appraisal data
+        housing_rec.load_remote_input_dataset("census_appraisal_data", census_appr_id)
+
+        housing_rec.set_parameter("result_name", result_name)
+
+        housing_rec.run_analysis()
+
+
+if __name__ == '__main__':
+    chained = False
+
+    run_with_base_class(chained)
diff --git a/tests/test_format.py b/tests/test_format.py
index a41de7708..84e302b0a 100644
--- a/tests/test_format.py
+++ b/tests/test_format.py
@@ -58,7 +58,10 @@
     os.path.join(PYINCORE_ROOT_FOLDER,
                  'tests/pyincore/analyses/joplinempiricalrestoration/test_joplinempiricalrestoration.py'),
     os.path.join(PYINCORE_ROOT_FOLDER, 'tests/pyincore/analyses/populationdislocation/test_populationdislocation.py'),
-    os.path.join(PYINCORE_ROOT_FOLDER, 'tests/pyincore/utils/test_popdisloutputprocess.py')
+    os.path.join(PYINCORE_ROOT_FOLDER, 'tests/pyincore/utils/test_popdisloutputprocess.py'),
+    os.path.join(PYINCORE_ROOT_FOLDER, 'pyincore/analyses/housingrecovery/housingrecovery.py'),
+    os.path.join(PYINCORE_ROOT_FOLDER, 'pyincore/analyses/housingrecovery/housingrecoveryutil.py'),
+    os.path.join(PYINCORE_ROOT_FOLDER, 'tests/pyincore/analyses/housingrecovery/test_housingrecovery.py')
 ]