Data Source

Data from 1994 - 2015
- U.S. GDP | Fred Economic Data | St.Louis Fed Click Here
- U.S. Income | U.S. Census Bureau, Current Population Survey Click Here

Theory

A probabilistic law focusing on uncertainty.
- An economic variable : \[\mbox{individual } i \mbox{ at time } t \mbox{ has } X_t(\omega^i).\]
- The change of $X_t(\omega^i)$ follows a probabilistic law $\mathbb{P}$: \[X_t(\omega^i) \sim \mathbb{P}(x_i,t).\]
- The (stationary) aggregate follows the same law $\mathbb{P}$: \[\sum_{i=1}^N X_t(\omega^i) \sim \mathbb{P}\left(\sum_{i=1}^N x^i,t\right).\]
Infinite divisible law (de Finetti, 1929)
- \[ Y_t = \sum_{i=1}^N X_t^i \overset{d}{=} X_t^1 + X_t^2 + \cdots + X_t^N\] \[ \mathbb{P}(y,t) = \mathbb{P}(x^{1},t) \star \mathbb{P}(x^{2},t) \star \cdots \star \mathbb{P}(x^{N},t) \] where $\star$ stands for the convolution operation.
Evolution of the economy \[\begin{array}{l} \mathbb{P}(x^{1})\quad\mbox{Adam's law}\\ \mathbb{P}(x^{1}),\mathbb{P}(x^{2})\quad\mbox{Adam's and Eve's laws}\\ \vdots \quad \quad \quad \vdots \quad \quad \ddots\\ \mathbb{P}(x^{1}),\mathbb{P}(x^{2}),\ldots,\mathbb{P}(x^{N},t)\quad\mbox{Current individuals' laws} \end{array}\]
Parameters of individual $\mathbb{P}$-law may be different from those in the aggregation.
- Aggregatable (Heterogenous) parameters
  - $\alpha_i$: Probability (unit of time) of a significant growth event for individual $i$
  - $c(x)$: Size of a significant growth.
- Non-aggregatable (Homogenous) parameters
  - $\beta$: Probability (unit of time) of a trivial growth event for everyone.
  - $1/\beta$: Size of a trivial growth.
$\mathbb{P}$-law is identifiable dynamically in terms of $(\alpha,\beta,c(x))$. \[ \frac{\partial}{\partial t}\mathbb{P}(x,t)=-c(x) \frac{\partial\mathbb{P}(x,t)}{\partial x}+\alpha\left[\int\left(\mathbb{W}(x|x')\mathbb{P}(x',t)\right)dx'-\mathbb{P}(x,t)\right] \] $\mathbb{W}(x|x')$ contains the parameter $\beta$.
- \[\mbox{Equilibrium solution: } \mathbb{P}(x)=\frac{\beta^{\alpha}x^{\alpha-1}e^{-\beta x}}{\Gamma(\alpha)}\sim\mbox{Gamma}(\alpha,\beta)\]
- \[\mbox{Aggregation: } X_{t}^{i} +X_{t}^{j} \sim\mbox{Gamma}(\alpha_{i}+\alpha_{j},\beta)\]
- \[\mbox{Evolution: } C\times \sum_{i=1}^{N} X_{t}^i\sim\mbox{Gamma}(\sum_{i=1}^{N}\alpha_i,\beta/C)\]
For the Gamma distribution: $\mathbb{E}[X] = \frac{\alpha}{\beta}$ and $\mbox{Var}[X] = \frac{\alpha}{\beta^2}$

$\alpha = (1,2,3,4,5),$ $\beta = (0.1,0.2,0.3,0.4,0.5),$ $\frac{\alpha}{\beta} = 10.$

$\alpha= 1,$ $\beta = (0.5,0.4,0.3,0.2,0.1),$ $\frac{\alpha}{\beta} = (2,2.5,3.3,5,10).$

$\alpha= (1,2,3,4,5),$ $\beta = 0.5,$ $\frac{\alpha}{\beta} = (2,4,6,8,10).$

$\alpha= (1,1.8,1.6,3.6,1.6),$ $\beta = (0.1,0.09,0.08,0.09,0.04).$ $\frac{\alpha}{\beta} = (10,20,20,40,40).$

Reduced Model

For a non-linear one-to-one function $g(\cdot)$ \[Y_{t}=g\left( X_{t}(\omega) \right),\] where $X_t(\omega)$ can represent the sum $\sum_{i=1}^{N}X_{t}(\omega^{i})$ or some $X_t(\omega)=(X_t(\omega^i),\dots,X_t(\omega^{j}))$ from individuals $i,\dots,j$.

$Y_{t}$ is the non-stationary aggregate: $Y_t \sim \mathbb{Q}(x,t).$
Mean field of $Y_t$: $m_{Y}(t)=\int y\mathbb{Q}(y,t)dy.$
Reduce form for this aggregation \[\begin{cases} \frac{dm_{Y}(t)}{dt} & = \mbox{Constant}+\mbox{Coef}_{1}m_{Y}(t)+\mbox{Coef}_{2}\sigma_{Y}^{2}(t) \\ \frac{d\sigma_{Y}^{2}(t)}{dt} & =\mbox{Constant}+\mbox{Coef}_{3}m_{Y}^{2}(t)+\mbox{Coef}_{4}\sigma_{Y}^{2}(t) \end{cases}\]
Growth 1994-2015: $m_Y(t)$

GDP = read_csv("~/Documents/2017/IncomeDistribution/data/A939RX0Q048SBEA.csv")
GDP = GDP[189:276,2]
GDP = ts(GDP,start=c(1994,1),frequency = 4)
plot(GDP, ylab="GDP = m_y(t)")

$\frac{dm_Y(t)}{dt}$

GDP.growth = diff(GDP)
GDP.growth[is.na(GDP.growth)]=0
plot(GDP.growth, ylab="d m_y(t)")

Construct $\sigma^2_Y(t)$ by AR residuals

GDP.ar = arima(GDP.growth, order = c(1,0,0))
GDP.ar

## 
## Call:
## arima(x = GDP.growth, order = c(1, 0, 0))
## 
## Coefficients:
##          ar1  intercept
##       0.3715   161.9089
## s.e.  0.0992    43.8295
## 
## sigma^2 estimated as 66911:  log likelihood = -606.86,  aic = 1219.71

hist(GDP.ar$resid, br=12)

sigma2 = GDP.ar$resid^2

GARCH(1,1) estimates the AR residual

GDP.garch = garch(GDP.ar$residuals, trace=FALSE)
summary(GDP.garch)

## 
## Call:
## garch(x = GDP.ar$residuals, trace = FALSE)
## 
## Model:
## GARCH(1,1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.4227 -0.4976  0.0188  0.6748  2.1834 
## 
## Coefficient(s):
##     Estimate  Std. Error  t value Pr(>|t|)  
## a0 6.092e+04   4.590e+04    1.327   0.1844  
## a1 1.748e-01   1.008e-01    1.734   0.0828 .
## b1 5.471e-10   5.981e-01    0.000   1.0000  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Diagnostic Tests:
##  Jarque Bera Test
## 
## data:  Residuals
## X-squared = 13.124, df = 2, p-value = 0.001413
## 
## 
##  Box-Ljung test
## 
## data:  Squared.Residuals
## X-squared = 1.6013, df = 1, p-value = 0.2057

plot(GDP.garch$residuals)

acf(na.omit(as.ts(GDP.garch$residuals)))

Mean field equaiton

\[ \frac{dm_{Y}(t)}{dt} = \mbox{Intercept}+\mbox{y[,2]}m_{Y}(t)+\mbox{y[,3]}\sigma_{Y}^{2}(t) \]

lagGDP.growth =lag(GDP.growth,1)
lagGDP.growth[is.na(lagGDP.growth)]=0
y = cbind(GDP.growth,GDP,sigma2)
GDP.lm = lm(y[,1]~ y[,2] + y[,3])
coef(summary(GDP.lm))

##                  Estimate   Std. Error   t value     Pr(>|t|)
## (Intercept)  6.651228e+02 3.067506e+02  2.168286 3.296532e-02
## y[, 2]      -9.628506e-03 6.705025e-03 -1.436014 1.547122e-01
## y[, 3]      -9.502302e-04 1.765819e-04 -5.381244 6.572483e-07

Residual analysis for the mean field equation

plot(GDP.lm)

ts.GDP.lm.r= ts(GDP.lm$residuals,start = c(1994), frequency = 4)
plot(ts.GDP.lm.r, col="gray", pch=16, xlab="", ylab="Mean field residuals")

qqnorm(ts.GDP.lm.r)
qqline(ts.GDP.lm.r,col=2)

acf(GDP.lm$residuals, main="")

Box.test(GDP.lm$residuals,1)

## 
##  Box-Pierce test
## 
## data:  GDP.lm$residuals
## X-squared = 0.24649, df = 1, p-value = 0.6196

garch.GDP.lm = garch(GDP.lm$residuals,trace = FALSE)
summary(garch.GDP.lm)

## 
## Call:
## garch(x = GDP.lm$residuals, trace = FALSE)
## 
## Model:
## GARCH(1,1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.2722 -0.5816 -0.0272  0.4895  3.7041 
## 
## Coefficient(s):
##     Estimate  Std. Error  t value Pr(>|t|)
## a0 4.996e+04          NA       NA       NA
## a1 1.006e-10          NA       NA       NA
## b1 6.293e-02          NA       NA       NA
## 
## Diagnostic Tests:
##  Jarque Bera Test
## 
## data:  Residuals
## X-squared = 17.798, df = 2, p-value = 0.0001365
## 
## 
##  Box-Ljung test
## 
## data:  Squared.Residuals
## X-squared = 0.034938, df = 1, p-value = 0.8517

Volatility of the mean field \[ \frac{d\sigma_{Y}^{2}(t)}{dt} =\mbox{Intercept}+\mbox{s[,2]} m_{Y}^{2}(t)+\mbox{s[,3]}\sigma_{Y}^{2}(t) \]

# Mean field volatility
sigma.dif = diff(sigma2)
sigma.dif[is.na(sigma.dif)]=0
s = cbind(sigma.dif, GDP^2, sigma2)
s.lm = lm(s[,1]~s[,2]+s[,3])
coef(summary(s.lm))

##                  Estimate   Std. Error     t value     Pr(>|t|)
## (Intercept) -7.306789e+03 9.991701e+04 -0.07312858 9.418797e-01
## s[, 2]      -2.121127e-05 4.686888e-05 -0.45256620 6.520419e-01
## s[, 3]       7.718935e-01 1.070106e-01  7.21324426 2.339573e-10

Residual analysis for the volatility equation

ts.s.lm.r= ts(s.lm$residuals,start = c(1994), frequency = 4)
plot(ts.s.lm.r, col="gray", pch=16, xlab="", ylab="Volatility residuals")

acf(s.lm$residuals, main="")

qqnorm(s.lm$residuals)
qqline(s.lm$residuals,col=2)

Box.test(s.lm$residuals,1)

## 
##  Box-Pierce test
## 
## data:  s.lm$residuals
## X-squared = 4.7398e-05, df = 1, p-value = 0.9945

new.garch = garch(GDP.lm$residuals,trace = FALSE)
summary(new.garch)

## 
## Call:
## garch(x = GDP.lm$residuals, trace = FALSE)
## 
## Model:
## GARCH(1,1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.2722 -0.5816 -0.0272  0.4895  3.7041 
## 
## Coefficient(s):
##     Estimate  Std. Error  t value Pr(>|t|)
## a0 4.996e+04          NA       NA       NA
## a1 1.006e-10          NA       NA       NA
## b1 6.293e-02          NA       NA       NA
## 
## Diagnostic Tests:
##  Jarque Bera Test
## 
## data:  Residuals
## X-squared = 17.798, df = 2, p-value = 0.0001365
## 
## 
##  Box-Ljung test
## 
## data:  Squared.Residuals
## X-squared = 0.034938, df = 1, p-value = 0.8517

Growth Rate (1994 - 2015): $\log m_Y(t) - \log m_Y(t)$

logGDP = log(GDP)
GDP.growth.rate = diff(logGDP) 
GDP.growth.rate[is.na(GDP.growth.rate)]=0
plot(GDP.growth.rate)

GDP.rate.ar = arima(GDP.growth.rate, order = c(1,0,0))
sigma2.rate = GDP.rate.ar$resid^2

# GARCH(1,1)

GDP.rate.garch = garch(GDP.rate.ar$residuals,trace=FALSE)
summary(GDP.rate.garch)

## 
## Call:
## garch(x = GDP.rate.ar$residuals, trace = FALSE)
## 
## Model:
## GARCH(1,1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -3.58454 -0.56414  0.02807  0.67152  2.41128 
## 
## Coefficient(s):
##     Estimate  Std. Error  t value Pr(>|t|)  
## a0 2.415e-05   1.597e-05    1.512   0.1304  
## a1 1.593e-01   8.016e-02    1.988   0.0468 *
## b1 3.397e-15   5.224e-01    0.000   1.0000  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Diagnostic Tests:
##  Jarque Bera Test
## 
## data:  Residuals
## X-squared = 8.0123, df = 2, p-value = 0.0182
## 
## 
##  Box-Ljung test
## 
## data:  Squared.Residuals
## X-squared = 0.95111, df = 1, p-value = 0.3294

plot(GDP.rate.garch$residuals)

acf(na.omit(as.ts(GDP.rate.garch$residuals)))

# Mean field equaiton

lagGDP.growth.rate =lag(GDP.growth.rate,1)
lagGDP.growth.rate[is.na(lagGDP.growth.rate)]=0
y = cbind(GDP.growth.rate,logGDP,sigma2.rate)
GDP.rate.lm = lm(y[,1]~ y[,2] + y[,3])
coef(summary(GDP.rate.lm))

##                 Estimate Std. Error   t value     Pr(>|t|)
## (Intercept)   0.17584136 0.07040236  2.497663 1.445226e-02
## y[, 2]       -0.01592906 0.00656432 -2.426613 1.737951e-02
## y[, 3]      -40.77699622 8.77836110 -4.645172 1.244946e-05

plot(GDP.rate.lm)

ts.GDP.rate.lm.r= ts(GDP.rate.lm$residuals,start = c(1994), frequency = 4)
plot(ts.GDP.rate.lm.r, col="gray", pch=16, xlab="", ylab="(log) Mean field residuals")

acf(GDP.rate.lm$residuals, main="")

Box.test(GDP.rate.lm$residuals, 1)

## 
##  Box-Pierce test
## 
## data:  GDP.rate.lm$residuals
## X-squared = 2.0689, df = 1, p-value = 0.1503

sigma.dif.rate = diff(sigma2.rate)
sigma.dif.rate[is.na(sigma.dif.rate)]=0
s = cbind(sigma.dif.rate, logGDP^2, sigma2.rate)
s.rate.lm = lm(s[,1]~s[,2]+s[,3])
coef(summary(s.rate.lm))

##                  Estimate   Std. Error     t value     Pr(>|t|)
## (Intercept)  4.050456e-05 4.436738e-04  0.09129357 9.274793e-01
## s[, 2]      -5.695058e-07 3.855121e-06 -0.14772709 8.829164e-01
## s[, 3]       7.800288e-01 1.071306e-01  7.28109908 1.721818e-10

ts.s.rate.lm.r= ts(s.rate.lm$residuals,start = c(1994), frequency = 4)
plot(ts.s.rate.lm.r, col="gray", pch=16, xlab="", ylab="(log) Volatility residuals")

acf(s.rate.lm$residuals, main="")

Box.test(s.rate.lm$residuals,1)

## 
##  Box-Pierce test
## 
## data:  s.rate.lm$residuals
## X-squared = 0.01114, df = 1, p-value = 0.9159

garch.rate = garch(GDP.rate.lm$residuals,trace=FALSE)
summary(garch.rate)

## 
## Call:
## garch(x = GDP.rate.lm$residuals, trace = FALSE)
## 
## Model:
## GARCH(1,1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -3.21105 -0.60261 -0.08058  0.46010  3.55524 
## 
## Coefficient(s):
##     Estimate  Std. Error  t value Pr(>|t|)
## a0 2.757e-05          NA       NA       NA
## a1 1.237e-02          NA       NA       NA
## b1 1.088e-09          NA       NA       NA
## 
## Diagnostic Tests:
##  Jarque Bera Test
## 
## data:  Residuals
## X-squared = 16.355, df = 2, p-value = 0.0002809
## 
## 
##  Box-Ljung test
## 
## data:  Squared.Residuals
## X-squared = 0.072048, df = 1, p-value = 0.7884

Structural Model

High order information can be exogenized by extracting the high order effects of $(Y_t,\dots,Y_t)$.
Structural form for this aggregation \[\begin{cases} m_{Y}(t) & =\mathbb{E}_{t}\left[g\left(\sum X_{t}^{i}\right)\mathbb{L}(t)\right]\\ \frac{dm_{Y}(t)}{dt} & =\theta_{t}\, m_{Y}(t)+e_{t} \\ X_{t} & \sim\mathbb{P}(x,t) \end{cases}\] where $\mathbb{L}(t)$ is the likelihood ratio for $Y_t$ and $X_t$; $\theta_{t}$ and $e_{t}$ are time varying constants.
Samples: People 15 years old and over
Measurement: Total money income
For information on confidentiality protection, sampling error, nonsampling error, and definitions, see www2.census.gov/programs-surveys/cps/techdocs/cpsmar16.pdf

DataUS <- read_delim("~/Documents/2017/IncomeDistribution/data/DataUS.txt", 
    "\t", escape_double = FALSE, trim_ws = TRUE)
names(DataUS)

##  [1] "Year Total"           "Total With Income"    "$1 to $2,499 or Loss"
##  [4] "$2,500 to $4,999"     "$5,000 to $7,499"     "$7,500 to $9,999"    
##  [7] "$10,000 to $12,499"   "$12,500 to $14,999"   "$15,000 to $17,499"  
## [10] "$17,500 to $19,999"   "$20,000 to $22,499"   "$22,500 to $24,999"  
## [13] "$25,000 to $27,499"   "$27,500 to $29,999"   "$30,000 to $32,499"  
## [16] "$32,500 to $34,999"   "$35,000 to $37,499"   "$37,500 to $39,999"  
## [19] "$40,000 to $42,499"   "$42,500 to $44,999"   "$45,000 to $47,499"  
## [22] "$47,500 to $49,999"   "$50,000 to $52,499"   "$52,500 to $54,999"  
## [25] "$55,000 to $57,499"   "$57,500 to $59,999"   "$60,000 to $62,499"  
## [28] "$62,500 to $64,999"   "$65,000 to $67,499"   "$67,500 to $69,999"  
## [31] "$70,000 to $72,499"   "$72,500 to $74,999"   "$75,000 to $77,499"  
## [34] "$77,500 to $79,999"   "$80,000 to $82,499"   "$82,500 to $84,999"  
## [37] "$85,000 to $87,499"   "$87,500 to $89,999"   "$90,000 to $92,499"  
## [40] "$92,500 to $94,999"   "$95,000 to $97,499"   "$97,500 to $99,999"  
## [43] "$100,000 and Over"

Income Data (1994 - 1996, 2013 - 2015): Numbers in thousands

us = unname(DataUS)
us = as.matrix(us)

head(us,3)

##      [,1]   [,2]  [,3]  [,4]  [,5]  [,6]  [,7]  [,8]  [,9] [,10] [,11]
## [1,] 1994 186402 18964 13353 17474 13683 14639 10750 11500  8656  9812
## [2,] 1995 188073 17946 12736 16638 13553 14389 10358 11855  8691 10383
## [3,] 1996 189997 16525 11984 16897 13358 14350 10551 11030  8740 10276
##      [,12] [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22]
## [1,]  6710  8106  5081  6748  3752  5283  2991  4487  2090  2625  1687
## [2,]  6864  8197  5247  7375  3976  5410  3092  4624  2264  2970  1923
## [3,]  7336  8784  5339  7664  4007  5677  3320  4898  2500  3091  2015
##      [,23] [,24] [,25] [,26] [,27] [,28] [,29] [,30] [,31] [,32] [,33]
## [1,]  2744  1336  1547  1022  1589   653   821   514   748   449   699
## [2,]  2871  1350  1813   990  1709   789  1100   586   894   417   832
## [3,]  3265  1485  1734   971  1956   889  1181   632   936   523   787
##      [,34] [,35] [,36] [,37] [,38] [,39] [,40] [,41] [,42] [,43]
## [1,]   357   542   214   298   298   208   338   179   195   128
## [2,]   376   601   347   333   180   356   207   212   138    72
## [3,]   444   689   313   479   288   369   237   266   166    97

tail(us,3)

##       [,1]   [,2]  [,3] [,4] [,5]  [,6]  [,7] [,8]  [,9] [,10] [,11] [,12]
## [20,] 2013 222003 14466 6764 9389 10931 12784 9529 11404  7966 10987  7346
## [21,] 2014 222972 14724 6613 8314 10795 13085 8898 10694  8101 11315  6894
## [22,] 2015 226762 14689 6262 7657 10551 12474 8995 10672  7931 11031  6962
##       [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23]
## [20,]  9047  5501  9462  4424  6983  4289  7824  3258  5116  3245  5939
## [21,]  9215  5640  9910  4234  7785  3941  8044  3042  5348  3113  6528
## [22,]  9623  5535 10399  4429  7975  3930  8091  3113  5718  3221  7130
##       [,24] [,25] [,26] [,27] [,28] [,29] [,30] [,31] [,32] [,33] [,34]
## [20,]  2249  4018  1984  4883  1559  2535  1316  3191  1215  2353  1051
## [21,]  2467  3968  1944  4806  1707  3022  1423  3227  1299  2753  1046
## [22,]  2489  3834  2066  5047  1894  3289  1493  3264  1372  2922  1307
##       [,35] [,36] [,37] [,38] [,39] [,40] [,41] [,42] [,43]
## [20,]  2213   982  1490   992  1742   725  1071   672 19108
## [21,]  2571  1001  1527   791  1738   683   976   723 19063
## [22,]  2725  1021  1508   856  1966   712  1090   768 20755

Ignore the year and truncate total number of individuals and income larger than $ 100,000

numcol = ncol(us)-1
numrow = nrow(us)
us = us[,3:numcol]

Estimate income distributions by the Gamma distribution

# Generate counting observations
alpha = matrix(0, numrow, 1)
beta = matrix(0,numrow,1)
sd.alpha = matrix(0, numrow, 1)
parm = matrix(0, numrow, 1)

t = 0

for (year in 1:numrow){
  t = t + 1
  iter = ncol(us)

  # For each year
  x.t = c()
  for (i in 1:iter){
    x = replicate(us[t,i], i)
    x.t = c(x.t,x) 
    }
  est.x = fitdist(x.t,"gamma") # MLE of gamma distribution
  parm[t] = est.x$estimate["shape"]/est.x$estimate["rate"]
  x.t = c()
  plot.legend = c("Year",1993+t)
  denscomp(list(est.x), xlab="", xlim=c(0,iter), fitcol = c("blue"), main = plot.legend, addlegend = FALSE)
  cdfcomp(list(est.x), xlab="", xlim=c(0,iter), fitcol = c("blue"), main = NULL, addlegend = FALSE);
  alpha[t] = est.x$estimate["shape"]
  beta[t] = est.x$estimate["rate"]
  sd.alpha[t] = est.x$sd
}

Dynamics of $\alpha_{t}$

up.bound = alpha + sd.alpha
low.bound = alpha - sd.alpha

plot(alpha~c(1994:2015),col="gray", pch=16, xlab="")
lines(up.bound, col=1, type="l", lty=3)
lines(low.bound, col=1, type="l", lty=3)

# abline(h=1.5, col=1)

Dynamics of $\beta_{t}$

plot(beta~c(1994:2015), col="gray", pch=16, xlab="")

Dynamics of $\frac{\alpha}{\beta}$ (mean)

plot(parm~c(1994:2015), col="gray", pch=16, xlab="", ylab="alpha/beta")

library(quantmod)
library(calibrate)

$(\alpha,\beta)$ give the endogenous forces to the growth \[\mathbb{E}\left[\sum_{i=1}^{N}X_{t}(\omega^{i})\right] =\frac{\alpha_{t}}{\beta_{t}}=\mbox{Cons}+\mbox{Coef}_{1}\times\ln m_{Y}(t)\]

GDP = na.omit(getSymbols("A939RX0Q048SBEA", src = "FRED", auto.assign = FALSE))
GDP = GDP['1993/2015']
GDP = ts(GDP,start=c(1993),frequency = 4)
GDP.sub = GDP[.indexmon(GDP)==9]
GDP.sub = ts(GDP.sub,start = c(1993), frequency = 1)
reg.pag = lm(parm~GDP.sub[-1])
plot(GDP.sub[-1],parm, col="gray", pch=16, ylab = "alpha/beta", xlab = "GDP per capita")
textxy(GDP.sub[-1],parm,c(1994:2015))
abline(coef=coef(reg.pag))

Endogeneity

\[ \hat{\varepsilon}_t = \mbox{Intercept} + \mbox{Coef} \times \frac{\alpha_t}{\beta_t} \]

reg.prim.endo = lm(reg.pag$resid~parm)
summary(reg.prim.endo)

## 
## Call:
## lm(formula = reg.pag$resid ~ parm)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.52909 -0.19375  0.03147  0.16826  0.48898 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.47040    0.51000  -0.922    0.367
## parm         0.04132    0.04450   0.928    0.364
## 
## Residual standard error: 0.2732 on 20 degrees of freedom
## Multiple R-squared:  0.04132,    Adjusted R-squared:  -0.006616 
## F-statistic: 0.862 on 1 and 20 DF,  p-value: 0.3643

$(\alpha,\beta)$ dynamics and the growth rate

GDP.growth.rate = diff(log(GDP.sub)) 
reg.pagr = lm(parm~GDP.growth.rate)
summary(reg.pagr)

## 
## Call:
## lm(formula = parm ~ GDP.growth.rate)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.4781 -0.7448  0.1578  0.8800  1.7620 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      11.9160     0.3477  34.272   <2e-16 ***
## GDP.growth.rate -35.6953    15.5002  -2.303   0.0322 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.22 on 20 degrees of freedom
## Multiple R-squared:  0.2096, Adjusted R-squared:  0.1701 
## F-statistic: 5.303 on 1 and 20 DF,  p-value: 0.03215

reg.endo = lm(reg.pagr$resid~parm)
plot(parm, reg.pagr$resid, col="gray",pch=16, ylab = "Residuals of Growth rate", xlab =  "alpha/beta" )
textxy(parm, reg.pagr$resid, c(1994:2015))
abline(coef=coef(reg.endo))

Endogeneity

\[ \hat{\varepsilon}_t = \mbox{Intercept} + \mbox{Coef} \times \frac{\alpha_t}{\beta_t} \]

reg.endo = lm(reg.pagr$resid~parm)
summary(reg.endo)

## 
## Call:
## lm(formula = reg.pagr$resid ~ parm)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.64393 -0.06395  0.15498  0.33365  0.63271 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -8.99869    1.04297  -8.628 3.56e-08 ***
## parm         0.79041    0.09101   8.685 3.20e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5587 on 20 degrees of freedom
## Multiple R-squared:  0.7904, Adjusted R-squared:  0.7799 
## F-statistic: 75.42 on 1 and 20 DF,  p-value: 3.204e-08

Estimate the growth rate function $\theta(t)$ \[\ln m_{Y}(t+1)-\ln m_{Y}(t) =\theta(t)+\varepsilon_{1,t}+\frac{\varepsilon_{2,t}}{\sum_{i=1}^{N}X_{t}(\omega^{i})}.\]
Filter the endgenous mean $\alpha/\beta$ and variance $\alpha/\beta^2$.

# Filter function
filtering = function(y,V,m0,C0){
  n  = length(y)
  a  = rep(0,n)
  R  = rep(0,n)
  m  = rep(0,n)
  C  = rep(0,n)
  for (t in 1:n){
    if(t==1){
      a[1] = m0
      R[1] = C0 
    }
    else{
      a[t] = m[t-1]
      R[t] = C[t-1] 
    }
    m[t]   = a[t]+(R[t]/(R[t] + V[t]))*(y[t]-a[t])
    C[t]   = R[t]*V[t]/(R[t] + V[t])
  }
  return(list(m=m,C=C))
}

# Estimation

m0 = GDP.growth.rate[1]
C0 = 1

t  = length(GDP.growth.rate)
theta  = matrix(0,t,1)
V = 1/(parm/beta)

estimate = filtering(GDP.growth.rate,V,m0,C0)    
theta = estimate$m
C = estimate$C
L  = theta - 1.96*sqrt(C)
U  = theta + 1.96*sqrt(C)

Endogeneity vanish.

resid.f   = GDP.growth.rate - theta
reg.res = lm(resid.f~parm)
summary(reg.res)

## 
## Call:
## lm(formula = resid.f ~ parm)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.044721 -0.005027  0.004239  0.008245  0.016058 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)
## (Intercept)  0.022986   0.027492   0.836    0.413
## parm        -0.002483   0.002399  -1.035    0.313
## 
## Residual standard error: 0.01473 on 20 degrees of freedom
## Multiple R-squared:  0.05084,    Adjusted R-squared:  0.00338 
## F-statistic: 1.071 on 1 and 20 DF,  p-value: 0.313

Filtered growth rate $\theta(t)$

plot(GDP.growth.rate[-1]~c(1995:2015), xlab ="", ylab="Growth Rate", col="gray",pch=16, ylim=c(ymin=-0.1, 0.1))
lines(theta~c(1994:2015),col=1)
lines(L~c(1994:2015),col=1, type="l", lty=3)
lines(U~c(1994:2015),col=1, type="l", lty=3)

Residual analysis

plot(parm, resid.f, xlab="alpha/beta", ylab="Residual of Structural Estimate", col="gray", pch=16)
textxy(parm, resid.f, c(1994:2015))
abline(coef=coef(reg.res))

Filtered Growth Rate (1994 - 2015)

Monte Carlo Updates (not included in the paper).

GDP.mc = GDP[.indexmon(GDP)==9]
GDP.mc = ts(GDP.mc,start = c(1993), frequency = 1)
GDP.gr.mc = diff(log(GDP.mc)) 

niter = 100000
V1    = V
for (iter in 1:niter){
    estimate.mc   = filtering(GDP.gr.mc,V1,m0,C0)
    theta.mc     = estimate.mc$m
    V1    = 1/rgamma(t,alpha+t/2,beta+sum((GDP.gr.mc-theta.mc)^2)/2)
}


estimate.mc = filtering(GDP.growth.rate,V1,m0,C0)    
theta.mc = estimate.mc$m
C.mc = estimate.mc$C
L.mc  = theta - 1.96*sqrt(C.mc)
U.mc  = theta + 1.96*sqrt(C.mc)

color = rgb(runif(1), runif(1), runif(1), .5)
plot(GDP.growth.rate[-1]~c(1995:2015), xlab="", ylab="Growth Rate (MC)", col="gray", pch=16, ylim=c(ymin=-0.1, 0.1))
lines(theta.mc~c(1994:2015),col=1)
lines(L.mc~c(1994:2015),col=1, type="l", lty=3)
lines(U.mc~c(1994:2015),col=1, type="l", lty=3)

resid.fmc = GDP.growth.rate - theta.mc
acf(resid.fmc, main="")

Box.test(resid.fmc,1)

## 
##  Box-Pierce test
## 
## data:  resid.fmc
## X-squared = 3.0917, df = 1, p-value = 0.07869

reg.res.mc = lm(resid.fmc~parm)
summary(reg.res.mc)

## 
## Call:
## lm(formula = resid.fmc ~ parm)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.044177 -0.005604  0.004082  0.008559  0.015211 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)
## (Intercept)  0.022215   0.027477   0.808    0.428
## parm        -0.002537   0.002398  -1.058    0.303
## 
## Residual standard error: 0.01472 on 20 degrees of freedom
## Multiple R-squared:  0.05302,    Adjusted R-squared:  0.005672 
## F-statistic:  1.12 on 1 and 20 DF,  p-value: 0.3026

plot(parm, resid.fmc, xlab="alpha/beta", ylab="Residual of Structural Estimate (MC)", col="gray", pch=16)
textxy(parm, resid.fmc, c(1994:2015))
abline(coef(reg.res.mc))

plot(c(1994:2015), theta.mc, xlab="", col="gray",pch=16, type="o")
abline(h=mean(theta))

Online Appendix: Uncertainty in Economic Growth and Inequality

Zhengyuan Gao (UCLouvain, CORE)

Data Source

Theory

Reduced Model

Structural Model