Most economic data are not generated by an experiment. They arrive with dependence across time, across space, and across the variables themselves. This dependence comes from market equilibrium, optimization, and selection. Each of those three forms of dependence violates the independence assumption of standard regression and complicates causal interpretation.
This chapter asks three questions.
Many observational datasets have temporal dependence, meaning that values at one point in time are related to past values. This violates the standard assumption of independence used in many statistical methods.
Stock prices are classic examples of temporally dependent processes. If Apple’s stock was high yesterday, it is more likely (but not guaranteed) to be high today.
# Load a CSV of daily Apple stock prices and plot the daily high as a line.
library('plotly')
stock <- read.csv('https://raw.githubusercontent.com/plotly/datasets/master/finance-charts-apple.csv')
fig <- plot_ly(stock, type = 'scatter', mode = 'lines')%>%
add_trace(x = ~ Date, y = ~ AAPL.High) %>%
layout(showlegend = FALSE)
fig# Today's high is not independent of past highs.A random walk is the simplest mathematical model of temporal dependence. Each new value is just the previous value plus a random shock (white noise).
# Generate Random Walk
tN <- 200
y_rw_sim <- numeric(tN)
y_rw_sim[1] <- stock$AAPL.High[1]
for (ti in 2:tN) {
y_rw_sim[ti] <- y_rw_sim[ti-1] + runif(1, -10, 10)
}
#x_sim <- runif(tN, -1, 1) White Noise
xy_rw_sim <- data.frame(Date=1:tN, RandomWalk=y_rw_sim)
fig <- plot_ly(xy_rw_sim, type = 'scatter', mode = 'lines') %>%
add_trace(x= ~ Date, y= ~ RandomWalk) %>%
layout(showlegend = FALSE)
fig# Today's value is not independent of past values.In a temporally dependent series, today’s value is not independent of past values. In contrast to cross-sectional data (e.g. individual incomes), time series often require special methods to account for memory and nonstationarity.
It often helps to see the marginal distribution alongside the time path. A nonstationary series can disguise itself in a histogram because the histogram pools across \(t\).
# Time path beside the marginal distribution
layout(matrix(c(1, 2), nrow = 1), widths = c(4, 1))
par(mar=c(4, 4, 1, 1))
plot(1:tN, y_rw_sim, type='l', lwd=2,
xlab='t', ylab=expression(y[t]))
title('Random walk: time path', font.main=1, adj=0)
par(mar=c(4, 1, 1, 2))
f_hat <- hist(y_rw_sim, breaks=20, plot=FALSE)
barplot(f_hat$counts, horiz=TRUE, space=0, border=NA,
col=grey(.5, .5), axes=FALSE)
axis(1)
title('Marginal', font.main=1, adj=0)
The simplest time-series methods implicitly assume the data have a stable distribution over time; this assumption fails often enough that we need a name for it.
Stationarity is useful as a precondition for most classical time-series tools. Note that there is more than one way to be nonstationary. A series with a deterministic trend is mean-nonstationary, while a random walk is variance-nonstationary because \(\mathbb{V}[Y_{t}]\) grows linearly in \(t\).
For example, consider the data generating process \(Y_{t} = \beta t + \epsilon_{t}\), where \(\epsilon_{t}\) is Normal with \(\mathbb{E}[\epsilon_{t}]=0\) and \(\sqrt{\mathbb{V}[\epsilon_{t}]}=\sigma+\alpha t>0\). Here \(\sigma\) is the baseline standard deviation in this Normal model. Then \(\mathbb{E}[Y_{t}]=\beta t\) and \(\mathbb{V}[Y_{t}]=(\sigma+\alpha t)^2\), so \(\beta\) affects mean stationarity and \(\alpha\) affects variance stationarity.
## MAIN EXAMPLE
tN <- 200
simulate_series <- function(beta, alpha, sigma=.2){
y_sim <- numeric(tN)
for (ti in 1:tN) {
mean_ti <- beta*ti
sd_ti <- (sigma + alpha*ti)
y_sim[ti] <- mean_ti + rnorm(1, sd=sd_ti)
}
return(y_sim)
}
# Plotting Functions
plot_setup <- function(alpha, beta){
plot.new()
plot.window(xlim=c(1, tN), ylim=c(-5, 20))
axis(1)
axis(2)
mtext(expression(y[t]), 2, line=2.5)
mtext('Time (t)', 1, line=2.5)
}
plot_title <- function(alpha, beta){
beta_name <- ifelse(beta==0, 'Mean Stationary', 'Mean Nonstationary')
alpha_name <- ifelse(alpha==0, 'Var Stationary', 'Var Nonstationary')
title(paste0(beta_name, ', ', alpha_name), font.main=1, adj=0)
}
par(mfrow = c(2, 2))
for(alpha in c(0, .015)){
for(beta in c(0, .05)){
plot_setup(alpha=alpha, beta=beta)
sim_cols <- c(rgb(.8, 0, 0, .5), rgb(0, 0, .8, .5))
for( sim in 1:2){
y_sim <- simulate_series(beta=beta, alpha=alpha)
lines(y_sim, col=sim_cols[sim], lwd=2)
}
plot_title(alpha=alpha, beta=beta)
}}
When values at one time depend on values at other times, we need correlation-based diagnostics across lags to characterize the dependence.
The ACF is useful for detecting temporal persistence (memory) within one series: for stationary processes it typically decays quickly, while nonstationary processes show ACFs that decay slowly or persist. The CCF is useful for detecting lagged relationships between two series, such as leading indicators or external drivers; if \(X_{t}\) is white noise, any visible structure in the CCF likely reflects nonstationarity in \(Y_{t}\).
## MAIN EXAMPLE
par(mfrow = c(2, 2))
for(alpha in c(0, .015)){
for(beta in c(0, .05)){
y_sim <- simulate_series(beta=beta, alpha=alpha)
acf(y_sim, main=NA)
plot_title(alpha=alpha, beta=beta)
}}
A cross-correlation plot diagnoses lagged dependence between two series at once:
## MAIN EXAMPLE
x_sim <- runif(tN, -1, 1) # White Noise
par(mfrow = c(2, 2))
for(alpha in c(0, .015)){
for(beta in c(0, .05)){
y_sim <- simulate_series(beta=beta, alpha=alpha)
ccf(y_sim, x_sim, main=NA)
plot_title(alpha=alpha, beta=beta)
}}
Many observational datasets exhibit spatial dependence, meaning that values at one location tend to be related to values at nearby locations. This violates the standard assumption of independent observations used in many classical statistical methods.
For example, elevation is spatially dependent: if one location is at high elevation, nearby locations are also likely (though not guaranteed) to be high. Similarly, socioeconomic outcomes like disease rates or income often cluster geographically due to shared environmental or social factors.
Just as stock prices today depend on yesterday, spatial variables often depend on neighboring regions, creating a need for specialized statistical methods that account for spatial autocorrelation.
Spatial data typically comes in two formats, each suited to different types of information:
# Vector Data
library('sf')
northcarolina_vector <- st_read(system.file('shape/nc.shp', package='sf'))
## Reading layer `nc' from data source `/home/jadamson/R-Libs/sf/shape/nc.shp' using driver `ESRI Shapefile'
## Simple feature collection with 100 features and 14 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
## Geodetic CRS: NAD27
plot(northcarolina_vector['BIR74'], main=NA)
title('Number of Live Births in 1974', font.main=1)
# https://r-spatial.github.io/spdep/articles/sids.html## MAIN EXAMPLE
# Raster Data
library('terra')
luxembourg_elevation_raster <- rast(system.file('ex/elev.tif', package='terra'))
plot(luxembourg_elevation_raster)
The size and shape of the areal units used to aggregate spatial data are a choice, and results can depend on it (Gehlke and Biehl 1934; Openshaw 1984; Fotheringham and Wong 1991), including for econometric estimates (Briant et al. 2010). See Adamson (2022) for an econometric model of ethnic maps and the endogenous areal unit problem, where the same variables shape both the units and the outcomes. Earlier work in geography studied units that are grouped non-randomly, but not units that form endogenously (Holt et al. 1996).
Just as with temporal data, stationarity in spatial data means that the statistical properties (like mean, variance, or spatial correlation) are roughly the same across space.
# Simulated 2D spatial fields
n <- 20
x <- y <- seq(0, 1, length.out = n)
grid <- expand.grid(x = x, y = y)
# 1. Stationary: Gaussian with constant mean and var
z_stationary <- matrix(rnorm(n^2, 0, 1), n, n)
# 2. Nonstationary: Mean increases with x and y
z_nonstationary <- outer(x, y, function(x, y) 3*x*y) + rnorm(n^2, 0, 1)
par(mfrow = c(1, 2))
# Stationary field
image(x, y, z_stationary,
main = NA,
col = terrain.colors(100),
xlab = 'x', ylab = 'y')
title('Stationary Field', font.main = 1)
# Nonstationary field
image(x, y, z_nonstationary,
main = NA,
col = terrain.colors(100),
xlab = 'x', ylab = 'y')
title('Nonstationary Field', font.main = 1)
Just as temporal data can exhibit autocorrelation, so can spatial data. Nearby locations tend to be more (or less) similar than spatial independence would predict.
Moran’s I is useful as the spatial analog of the ACF. It summarizes overall spatial dependence in one number. Large positive values signal that nearby observations share information, so we cannot treat them as independent. A two-variable analog measures spatial cross-correlation, computed like the CCF but with spatial rather than temporal lags.
## MAIN EXAMPLE
# Raster Data Example
autocor(luxembourg_elevation_raster, method='moran', global=TRUE)
## elevation
## 0.8917057Cross-Correlation. We can also assess the relationship between two variables at varying distances.
# Vector Data Example
xy_hat <- as.data.frame(northcarolina_vector)[, c('BIR74', 'SID74')]
column_means <- colMeans(xy_hat)
# Format Distances
dmat <- st_distance( st_centroid(northcarolina_vector) )
dmat <- units::set_units(dmat, 'km')
# At Which Distances to Compute CCF
# summary(dmat[, 1])
rdists <- c(-1, seq(0, 100, by=25)) # includes 0
rdists <- units::set_units(rdists , 'km')
# Compute Cross-Covariances
varXY <- prod( apply(xy_hat, 2, sd) )
CCF <- lapply( seq(2, length(rdists)), function(ri){
# Which Observations are within (rmin, rmax] distance
dmat_r <- dmat
d_id <- (dmat_r > rdists[ri-1] & dmat_r <= rdists[ri])
dmat_r[!d_id] <- NA
# Compute All Covariances (Stationary)
covs_r <- lapply(1:nrow(dmat_r), function(i){
pairsi <- which(!is.na(dmat_r[i, ]))
covXiYj <- sapply(pairsi, function(j) {
dXi <- xy_hat[i, 1] - column_means[1]
dYj <- xy_hat[j, 2] - column_means[2]
return(dXi*dYj)
})
return(covXiYj)
})
corXY <- unlist(covs_r)/varXY
return(corXY)
} )# Plot Cross-Covariance Function
x <- as.numeric(rdists[-1])
par(mfrow=c(1, 2))
# Distributional Summary
# The axis stops just past -1 and 1, so 106 of the 1708 values are not drawn
boxplot(CCF,
border=NA, medcol=grey(.4), # a darker shade of the default grey fill, for the median
outline=FALSE, whisklty=0, staplelty=0,
ylim=c(-1, 1), #quantile(unlist(CCF), probs=c(.05, .95)),
names=x,
main=NA,
xlab='Distance [km]',
ylab='Cross-Correlation of BIR74 and SID74')
# Raw Observations, in place of whiskers and outliers
stripchart(CCF,
pch='-', col=grey(0, .25), cex=1,
vert=TRUE, add=TRUE)
title('Binned Medians and IQRs', font.main=1, adj=0)
abline(h=0, col=grey(.5), lty=3)
# Inferential Summary
CCF_means <- sapply(CCF, mean)
plot(x, CCF_means,
ylim=c(-1, 1),
type='o', pch=16,
main=NA,
xlab='Distance [km]',
ylab='Cross-Correlation of BIR74 and SID74')
title('Binned Means + 95% Band', font.main=1, adj=0)
abline(h=0, col=grey(.5), lty=3)
# Quick and Dirty Subsampling CI
CCF_meanCI <- sapply(CCF, function(corXY){
ss_size <- floor(length(corXY)*3/4)
corXY_boot <- sapply(1:200, function(b){
corXY_b <- sample(corXY, ss_size, replace=FALSE)
mean(corXY_b, na.rm=TRUE)
})
quantile(corXY_boot, probs=c(.025, .975), na.rm=TRUE)
})
polygon( c(x, rev(x)),
c(CCF_meanCI[1, ], rev(CCF_meanCI[2, ])),
col=grey(0, .25), border=NA)
Two variables can also be related locally, within a window around each location (Brunsdon et al. 2002; Adamson 2021, 2025). Neighbors’ values can also enter a regression directly: a spatially lagged X model regresses outcomes on each location’s own variables and the sum of its neighbors’ values (Anselin 1988; LeSage and Pace 2009; Halleck Vega and Elhorst 2015; Adamson 2025).
Spatial dependence also matters for standard errors (Conley 1999; Müller and Watson 2022; Conley and Kelly 2025). For example, Adamson (2025) uses a spatial subsampling procedure (Politis et al. 1999; Sherman 1996; Choi and Chue 2007) for strongly dependent data.
In addition to spatial and temporal dependence, many observational datasets exhibit interdependence between variables for economic reasons.
Economic outcomes are not free variables we can move at will. Many are jointly determined, affected by selection, or driven by variables we have not measured.
Endogeneity is a common reason that observational data are hard to interpret causally. In the linear model \(Y_{i}=\sum_{k=0}^{K}\beta_{k}X_{ik} + \epsilon_{i}\), endogeneity means \(\mathbb{E}[X_{ik}\epsilon_{i}] \neq 0\) for some \(k\). OLS is then biased away from the structural causal parameter.1
## MAIN EXAMPLE
# Simulate data with an endogeneity issue
n <- 300
z_sim <- rbinom(n, 1, .5)
xy_sim <- sapply(z_sim, function(zi){
y_sim <- rnorm(1, zi, 1)
x_sim <- rnorm(1, zi*2, 1)
c(x_sim, y_sim)
})
xy_sim <- data.frame(x_sim=xy_sim[1, ], y_sim=xy_sim[2, ])
plot(y_sim ~ x_sim, data=xy_sim, pch=16, col=grey(0, .25), main=NA)
abline(lm(y_sim ~ x_sim, data=xy_sim), col=rgb(1, 0, 0, .8), lwd=2)
I will focus on the seminal economic example to provide some intuition.
This model has three structural relationships. First, market supply is the sum of quantities supplied by individual firms at a given price. Second, market demand is the sum of quantities demanded by individual people at a given price. Third, market supply equals market demand in equilibrium. Assuming market supply and demand are linear, we can write these three relationships as \[\begin{aligned} Q_{S}(P) &= \alpha_{S} + \beta_{S} P + E_{S}, \\ Q_{D}(P) &= \alpha_{D} - \beta_{D} P + E_{D}, \\ Q_{D} &= Q_{S} = Q. %% $Q_{D}(P) = \sum_{i} q_{D,i}(P)$, \end{aligned}\]
Neither variable is “exogenous” when the two are jointly determined in equilibrium, since we can only observe their joint outcomes.
The reduced form is useful for seeing what observational data can identify. The reduced-form intercepts and shocks combine the structural parameters in known ways. Recovering the structural slopes \(\beta_{S}, \beta_{D}\) requires additional identifying variation. The equilibrium condition \(Q_{D} = Q_{S}\) in the structural equations above makes price and quantity jointly determined, so both are outcomes.
## MAIN EXAMPLE
# Demand Curve Simulator
qd_fun <- function(p, Ad=8, Bd=-.8, Ed_sigma=.25){
Qd <- Ad + Bd*p + rnorm(1, 0, Ed_sigma)
return(Qd)
}
# Supply Curve Simulator
qs_fun <- function(p, As=-8, Bs=1, Es_sigma=.25){
Qs <- As + Bs*p + rnorm(1, 0, Es_sigma)
return(Qs)
}
# Quantity Supplied and Demanded at 3 Prices
cbind(P=8:10, D=qd_fun(8:10), S=qs_fun(8:10))
## P D S
## [1,] 8 1.8129598 0.1162905
## [2,] 9 1.0129598 1.1162905
## [3,] 10 0.2129598 2.1162905
# Market Equilibrium Finder
eq_fun <- function(demand, supply, P){
# Compute EQ (what we observe)
eq_id <- which.min( abs(demand-supply) )
eq <- c(P=P[eq_id], Q=demand[eq_id])
return(eq)
}## MAIN EXAMPLE
# Simulations Parameters
N <- 300 # Number of Market Interactions
P <- seq(5, 10, by=.01) # Price Range to Consider
# Generate Data from Competitive Market
# Plot Underlying Process
plot.new()
plot.window(xlim=c(0, 2), ylim=range(P))
EQ1 <- sapply(1:N, function(n){
# Market Data Generating Process
demand <- qd_fun(P)
supply <- qs_fun(P)
eq <- eq_fun(demand, supply, P)
# Plot Theoretical Supply and Demand
lines(demand, P, col=grey(0, .01))
lines(supply, P, col=grey(0, .01))
points(eq[2], eq[1], col=grey(0, .05), pch=16)
# Save Data
return(eq)
})
axis(1)
axis(2)
mtext('Quantity', 1, line=2)
mtext('Price', 2, line=2)
Suppose we ask “what is the effect of price on quantity?” You can simply run a regression of quantity, \(\hat{y}\), on price, \(\hat{x}\), and get the slope coefficient \(\hat{b}_{1} = \hat{c}_{Q^{*} P^{*}} / \hat{v}_{P^{*}}\). You always get a number back, but it is hard to interpret meaningfully.
## MAIN EXAMPLE
# Analyze Market Data
xy1_sim <- data.frame(t(EQ1), cost='1', T=1:N)
reg1 <- lm(Q ~ P, data=xy1_sim)
summary(reg1)
##
## Call:
## lm(formula = Q ~ P, data = xy1_sim)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.47571 -0.12130 -0.00619 0.11426 0.48781
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.48931 0.43247 1.131 0.259
## P 0.04315 0.04867 0.887 0.376
##
## Residual standard error: 0.1715 on 298 degrees of freedom
## Multiple R-squared: 0.002631, Adjusted R-squared: -0.0007155
## F-statistic: 0.7862 on 1 and 298 DF, p-value: 0.376This simple example has a profound insight: price-quantity data does not generally tell you how price affects quantity (or vice-versa). The reason is simultaneity: price and quantity mutually cause one another in markets.2
Moreover, this example also clarifies that our initial question “what is the effect of price on quantity?” is misguided. We could more sensibly ask “what is the effect of price on quantity supplied?” or “what is the effect of price on quantity demanded?”
With multiple linear regression, endogeneity biases are not just a problem for your main variable of interest. Suppose you are interested in how \(X_{i1}\) affects \(Y_{i}\), conditioning on \(X_{i2}\), and that the data generating process is linear: \(Y_{i}=\beta_{0}+\beta_{1}X_{i1}+\beta_{2}X_{i2}+\epsilon_{i}\). You paid special attention in your research design to find a case where \(X_{i1}\) is truly exogenous. Unfortunately, if \(X_{i2}\) is correlated with the error term, then there is also a bias for \(X_{i1}\). The magnitude of the bias for \(X_{i1}\) depends on the correlations between \(X_{i1}\) and \(X_{i2}\) as well as \(X_{i2}\) and \(\epsilon_{i}\).3
Comment the script you wrote for this chapter, then restart R and check that it runs from a clean session, then check the script with AI as explained in Working with AI. Write three sentences from memory on the main statistical idea of this chapter, and ask the assistant what is wrong, vague, or missing. Finish with your own questions about whatever you found hardest.
Explain why simultaneity in a competitive market makes it difficult to interpret the OLS slope \(\hat{b}_{1}\) from a regression of quantity on price. In your answer, distinguish between the structural demand slope \(\beta_{D}\) and what OLS actually estimates.
Using the North Carolina births data (sf::st_read(system.file('shape/nc.shp', package='sf'))), compute Moran’s I for the variable SID74 (sudden infant deaths in 1974). Is there evidence of spatial autocorrelation? Briefly explain what a high value of Moran’s I would imply for standard OLS inference.
Simulate a random walk of length 200 in R (starting at 0, with increments drawn from runif(1, -1, 1)). Plot the series, compute the ACF using acf(), and describe how the ACF pattern differs from what you would expect for white noise.
This chapter introduced temporal, spatial, and economic dependence in observational data. It measured temporal dependence with the ACF and CCF, and spatial dependence with Moran’s I. For economic dependence, it listed the sources of endogeneity. The supply-demand simulation made the endogeneity bite concrete: we ran 300 market periods through qd_fun/qs_fun, fit lm(Q ~ P, data=xy1_sim), and got a slope that confuses the supply and demand sides because price and quantity are jointly determined.
\(X_{ik}\) and \(\epsilon_{i}\) may be correlated for other reasons too, such as when \(X_{ik}\) is measured with error.↩︎
Although there are many ways this simultaneity can happen, economic theorists have made great strides in analyzing the simultaneity problem as it arises from equilibrium market relationships. In fact, 2SLS arose to understand agricultural markets. With a linear structure to supply and demand, we can even use algebra to solve for the equilibrium price and quantity analytically as \[\begin{aligned} P^{*} &= \frac{\alpha_{D}-\alpha_{S}}{\beta_{D}+\beta_{S}} + \frac{E_{D} - E_{S}}{\beta_{D}+\beta_{S}}, \\ Q^{*} &= \frac{\alpha_{S}\beta_{D}+ \alpha_{D}\beta_{S}}{\beta_{D}+\beta_{S}} + \frac{E_{S}\beta_{D}+ E_{D}\beta_{S}}{\beta_{D}+\beta_{S}}. \end{aligned}\]↩︎
Stack the rows \([X_{i1}, X_{i2}]\) into the matrix \(\mathbf{X}\). Stack the \(Y_{i}\) into the vector \(\mathbf{Y}\), and let \(\mathbf{B}\) be the OLS coefficient estimator. After algebraic work, you can find the bias \(\mathbb{E}[ \mathbf{B} - \beta]\) suffers from a contamination effect. \[\begin{aligned} \mathbf{B} &= [\mathbf{X}'\mathbf{X}]^{-1}\mathbf{X}'\mathbf{Y} \\ \mathbb{E}[\mathbf{X}'\epsilon] &= \begin{bmatrix} 0 \\ \rho \end{bmatrix} \\ \mathbb{E}[ \mathbf{B} - \beta] &= [\mathbf{X}'\mathbf{X}]^{-1} \begin{bmatrix} 0 \\ \rho \end{bmatrix} \neq \begin{bmatrix} 0 \\ \rho \end{bmatrix} \end{aligned}\] In words: when \(X_{i2}\)’s row of \(\mathbf{X}'\epsilon\) has expectation \(\rho\), the bias spreads through the inverse \((\mathbf{X}'\mathbf{X})^{-1}\) onto every coefficient. That includes the one you cared about, on \(X_{i1}\).↩︎