The Regime Shift Associated with the 2004–2008 US Housing Market Bubble

The Subprime Bubble preceding the Subprime Crisis of 2008 was fueled by risky lending practices, manifesting in the form of a large abrupt increase in the proportion of subprime mortgages issued in the US. This event also coincided with critical slowing down signals associated with instability, which served as evidence of a regime shift or phase transition in the US housing market. Here, we show that the US housing market underwent a regime shift between alternate stable states consistent with the observed critical slowing down signals. We modeled this regime shift on a universal transition path and validated the model by estimating when the bubble burst. Additionally, this model reveals loose monetary policy to be a plausible cause of the phase transition, implying that the bubble might have been deflatable by a timely tightening of monetary policy.


S1 Data
The data for the three US housing market variables were downloaded from Zillow (http:// www.zillow.com/research/data/), which is an online real estate database company in the US. Because of its popularity (Zillow claims that it has data on more than 110 million US homes), its data represents a considerable cross-section of the US housing market. The three variables that we downloaded from Zillow are the median sale price (transaction price, Jul 97 -Dec 12), homes sold (transaction volume, Jul 97 - Dec 12), and % of homes sold for gain (Jan 98 - Feb 13). The data for these variables come in the form of monthly time series for a number of cities in the US. When computing the average normalized trajectory, in order for a comparison across variables, we only use data from time periods which the three variables overlap i.e. Jan 98 -Dec 12. It should be noted that we downloaded the data in early 2013. Due to an update by Zillow to the methodology of computing some variables in 2014, some of the data we used might not be available on the website now. In this case, requests for data should be directed either to Zillow or the corresponding author using the e-mail provided on the main text.

S2 Early Warning Signals
We calculate the lag-1 autocorrelation, standard deviation, spectral reddening, and skewness for the homes sold variable in overlapping 24-month time windows that we slide one month at a time. A time window length of 24 months was chosen as a balance between sensitivity to events in the housing market and resolution of the statistical early warning signals. The detection of statistically significant trends in aggregated early warning signals is robust to the choice of time window length (Section S8). Only cities that have no missing data from the homes sold variable are used to calculate these early warning signals.
From the theory of nonlinear dynamical systems, decreasing stability in a stable attractor means that recovery from perturbations take longer, giving rise to an increasing autocorrelation, spectral reddening of the power spectrum and an increasing standard deviation [8,30]. Also, as the attractor becomes less stable, the energy well associated with it may become asymmetric [17]. In particular, the skewness also gives the direction in which the free energy is flattening. Let {r t }, with indices indicating the time index, represent a detrended time series window for the transaction volume variable of the normalized average trajectory for any particular city. Then the lag-1 autocorrelation is calculated by the formula R = E[(r t − µ)(r t+1 − µ)], where µ is the mean and σ 2 = E[(r t − µ) 2 ]/23 is the variance. Before calculating the standard deviation of a time window, the parent time series of that time window is normalized by dividing by the standard deviation of the parent time series. Spectral reddening is quantified by calculating the median frequency ω 1/2 of the frequency spectrum; a lower ω 1/2 indicates a frequency spectrum that is more redshifted, having more power in the lower frequencies. The frequency spectrum is computed using a discrete Fourier transform (DFT), with the constant discarded and the frequency spectrum normalized. Finally, the skewness is computed according to the formula Early warning signals for each time window were calculated. An early warning signal at any particular time was aggregated across all cities by counting the number of cities that pass a certain threshold for that early warning signal. This results in an aggregated early warning signal for the US housing market. The threshold was chosen such that the aggregated early warning signal was sensitive enough to the occurrence of statistically significant trends without the threshold being too stringent such that the trends were undetectable ( Figure A). The thresholds for the autocorrelation, standard deviation, spectral reddening and skewness are 0.5, 1.5, the third lowest frequency and 1.35 respectively.
Before calculating the statistical early warning signals, we detrend each time series using LOESS smoothing, after which individual time windows can be segmented from the time series of residues. This was done using the 'smooth' function in MATLAB with the 'rloess' method and a span of 40%. A span of 40% was chosen so that the smoothing fitted the slow dynamics well enough without overfitting the fast dynamics of the time series ( Figure B). The detection of statistically significant aggregated early warning signals around the period of the two transitions is robust to the choice of span (Section S8).  Significance Testing Significance testing of early warning signals is done by calculating the significance of the Kendall's τ [30], a measure of rank correlation between two random variables. In this case, one of the random variables is an aggregated early warning signal while the other is time. Hence, Kendall's τ measures the persistence of an increasing trend in early warning signal as the stability of the attractor decreases towards a phase transition point. The null model time series is created by bootstrapping from the time series of each city before the start of an increasing trend of aggregated early warning signal. For the spectral reddening and autocorrelation early warning signals, moving block bootstraps were used instead to retain the memory from the original time series [31]. We then calculated the early warning signal for the null model time series of each city and obtained the aggregated early warning signal which is the number of cities that passed a certain early warning threshold as described previously. After creating 10,000 such instances of aggregated early warning signal from null model time series, we can calculate the probability of obtaining a Kendall τ greater than or equal to that observed in the increasing trend. We consider the trend to be significant when the probability of obtaining such a Kendall τ due to random chance (null model) is less than 0.05. In Section S8, we conduct a sensitivity analysis on the span and the time window length and show that the detection of statistically significant aggregated early warning signals around the time of both transitions is robust.

S3 Average Normalized Trajectory
To study the trajectory of the US housing market in the phase space of the three variables, we partitioned each variable's time series into two-year sliding time windows for the purpose of a later comparison of the trajectory with the early warning signals. Each time window is normalized by dividing by the average of the entire time series, of which the time window is a constituent. The average normalized trajectory of the housing market for a variable at the time index of the time window is then calculated by averaging across the cities, the average of each city's time window for that variable. Doing this across all three variables and all time windows produces a time series which we used as the average normalized trajectory of the US housing market. Only cities that have no missing data from all three variables are used to calculate the average normalized trajectory.

S4 Fit Results of the Homes Sold Trajectory
The dynamical equation used in the main text is where U(x, t) is the free energy approximated by a Taylor expansion Before fitting, it is possible to remove the third-order order term C 1 by a fixed translation in x so that C 1 = 0. Hence, C(t) may be set to zero after the driving phase. We relate the homes sold trajectory Q(t) to x(t) by where  Determining Q g , the natural growth rate Because US population growth was roughly linear for the past two decades, we assumed that the natural growth in Q(t) is linear in time i.e. Q g t (Eq (D)). To estimate Q g and Q 0 , we perform a linear regression on Q(9 ≤ t ≤ t 0 ) which is the linear portion of Q(t) before the kink at t 0 when the driving was started. Here t is given in index units so t = 9 is the 9th time window in sequence. The linear regression gives the estimate  Fitting algorithm All fits were done with MATLAB's lsqnonlin function, a nonlinear regression routine which uses a trust-region algorithm. Numerical integration of Eq (B) was done using a fourth-order Runge-Kutta method with a step size of 10 −4 month. The parameters fitted were Sensitivity analysis To ensure our model is robust, we explored alternatives to the assumptions used in the model, namely, concave and convex functions for Q ∆ (t) instead of the linear ones used in the model. We also accounted for different values of Q g within the confidence interval estimated earlier. After fitting, our results Section S7 show that any fits with an R 2 larger than that of the original fit did not show any significant deviation of t ∆ that would affect the validity of the model used.
Endogeneity We used a linear regression for the time series preceding the transition so as to obtain a natural growth rate for the transaction volume in the United States. We assumed that this growth is linear because the population growth in the United States was roughly linear for the past few decades. In the Supporting Information Figure S3, we have plotted a regression and also a plot of residuals. Judging from the plot of residuals, we feel that the linear regression is satisfactory. Because we did sensitivity analysis for the parameters estimated from the linear regression, we feel that the model is robust even if there exists some explanatory variables that were omitted from the regression. For the regime shift model, it is safe to assume that not all endogenous variables are accounted for. However, since a large amount of variation in the homes sold variable can be explained by the model, we think that it is safer not to include more variables and run the risk of overfitting the data with variables that might only be spuriously related to the homes sold variable.  This rules out the possibility that Q(t) may be explained by an absence of driving in the model. The black curve shows an extrapolation of the fitted dynamics under a time reversal from the red curve. The local minimum in the black curve shows a little ghosting, a remnant of the initial fixed point that was annihilated by the saddle node bifurcation during the driving phase, slowing trajectories moving through the region of phase space near the former position of the initial fixed point.

S5 Fitting U (x, t) to simulated data of a grazing model
A commonly used one-dimensional grazing model in ecology exhibiting a regime shift is given by the dimensionless dynamical equation [32,33] Here, X is the vegetation biomass as a proportion of its carrying capacity, γ controls the grazing rate, α determines at what level of X the grazing rate saturates and τ is the dimensionless time.  (Table D) as an example, we demonstrate a soft landing without a phase transition by stopping the driving before the bifurcation occurs.
The model exhibits two stable fixed points when γ = 0.23 and α = 0.1. Increasing γ causes a saddle-node bifurcation at γ = 0.26 which eliminates the stable fixed point with the larger vegetation biomass. If the system is initially at the stable fixed point with the larger vegetation biomass, the system will undergo a regime shift to the remaining stable fixed point which is at the smaller vegetation biomass. We generate simulated data for X(τ ) by numerically integrating Eq F. α is kept at a value of 0.1. From 0 ≤ τ ≤ 10, the system is initially at the high biomass stable fixed point. When 10 < τ ≤ 40, γ(τ ) is increased linearly with τ from γ(10) = 0.23 to γ(40) = 0.28. Then for 40 < τ ≤ 80, γ is kept constant at 0.28. The simulated data X(40 < τ ≤ 80) was then fitted to Eq B with free energy given by Eq C. Here, A(40 < τ ≤ 80) = A 1 , B(40 < τ ≤ 80) = B 1 , C(40 < τ ≤ 80) = C 1 = 0 and D(40 < τ ≤ 80) = D 1 are kept constant for the fit. The fitted parameters are A 1 , B 1 , D 1 and x 1 , where x 1 is the initial position of the system at τ = 40.
The approximation of Eq B to the grazing model given by Eq F depends on how close the system is to the target regime. We modified the grazing model to improve this approximation. This was done by adding a new parameter k to the grazing model resulting in the modified model The parameters used in the modified model are k = 0.9 and α = 0.16. From 0 ≤ τ ≤ 10, the system is initially at the high biomass stable fixed point. When 10 < τ ≤ 80, γ(τ ) is increased linearly with τ from γ(10) = 0.231 to γ(80) = 0.236. Then for 80 < τ ≤ 160, γ(τ ) is kept constant at 0.236. This modification allows us to decrease the distance between the initial regime and the final regime without substantially increasing the coefficients of the 5th order and higher terms in Eq C thus improving the approximation. The results of the fit to simulated data from the original and modified model are shown in Table C. In Figure Ed, we contrast the behavior of the grazing model under extended forcing (resulting in biomass collapse), and when forcing is terminated early (allowing the biomass to achieve a soft landing).

S6 Overfitting the Driving Phase
In the main text, we alluded to the problem of overfitting the driving phase when estimating the functional forms of A(τ ), B(τ ), C(τ ) and D(τ ) without additional mechanistic insight. The problem lies in the large number of parameters involved relative to the simplicity of the curve being fitted. Additionally, error propagation from the estimation of A 1 , B 1 and D 1 after the driving phase is also a problem when trying to model the driving phase. As an example, consider the modified grazing model specified in Section S5 and Eq G. We first try to estimate A(τ ), B(τ ), C(τ ) and D(τ ) during the driving phase (10 < τ ≤ 80) by assuming that A(τ ), B(τ ), C(τ ) and D(τ ) are linear functions, where Y represents A, B, C or D. In fact, A(τ ), B(τ ), C(τ ) and D(τ ) are indeed linear functions during the driving phase in this example because γ(τ ) is driven linearly when generating the simulated data. We fit the driving phase of the simulated data to Eqs B, C and H after obtaining estimates of A 1 , B 1 and D 1 from fitting the portion of data after the driving phase (Table C). The parameters to be fitted are thus A 0 , B 0 , C 0 and D 0 . As we can see, very different parameters (Table D) fit the data almost equally well (Fig. F), and therefore we have a problem of the model overfitting the data.

S7 Sensitivity Analysis (Fitted Model)
The 95% confidence interval in the estimation of Q g =5.564E-4 is (3.654E-4, 7.743E-4). We conduct the sensitivity analysis with five different values of Q g : 3.654E-4, 4.609E-4, 5.564E-4, 6.518E-4 and 7.743E-4. We also explore different convex and concave functions for Q ∆ (t). For the convex functions, we use an exponential function of the form A comparison of the exact theoretical approximation (Eq B and C) with parameters calculated directly from the modified grazing model (Eq G) against two alternate local solutions obtained by fitting the approximation (Eq B and C) and the assumed linear driving forces (Eq H) to simulated data from the modified grazing model during the driving phase. A 1 , B 1 and C 1 were calculated from the modified grazing model for all three plots. The fits are both better than the theoretical approximation as far as R 2 is concerned, but their parameters are very different, and also very different from the theoretical approximation (Table D). Table D: Parameters recovered from two alternate local solutions to the fitting algorithm in addition to the R 2 and τ b , which is the time the bifurcation occurred during the driving phase.
For the convex exponential function, we fit for values of λ starting from 5E-3, increasing in increments of 5E-3 until 7E-2. For the convex power function, we fit for values of λ starting from 1.05, increasing in increments of 0.05 until 1.70. For the concave functions, we use the power function given by Eq J and the logarithmic function found by inverting Eq I and applying the appropriate linear transformations. For the concave power function, we fit for values of λ starting from 0.3, increasing in increments of 0.05 to 0.95. For the logarithmic function, we fit for values of λ starting from 5E-3, increasing in increments of 5E-3 until 7E-2. Hence, increasing λ is associated with an increasing curvature in Q ∆ (t) for the exponential, logarithmic and convex power function whereas increasing λ is associated with a decreasing curvature in Q ∆ (t) for the concave power function. The results of the sensitivity analysis may be found in Figure G. The distribution of t ∆ and its 95% confidence interval bounds for fits with an R 2 greater than the original fit can be found in Figure H. The results of the sensitivity analysis show that the estimation of t ∆ is robust for the alternative scenarios considered in the sensitivity analysis.

S8 Sensitivity Analysis (Early Warning Signals)
The following tables contain the sensitivity analysis of the p-values of increasing trends in aggregated early warning signals around the Subprime Loans Transition and the Financial Crisis Transition with respect to the span and time window lengths. The span is given as a percentage of the total number of data points in the time series. 1,000 samples were obtained to compute each p value except for the parameters used for presenting the results of the main text where the span is 40% and the time window length is 24. 10,000 samples were used in this case.