Panel Data: Cigarette
Alex Verville
2026-07-24
library(AER)
packages <- c("AER", # load dataset
"plm", # for FE, does what ivreg does with more controls
"stargazer", # summary stats for sharing
"ggplot2", #for graphing,
"psych",
"naniar", # for visualisation of missing data,
"visdat", # for visualisation of missing data,
"VIM", # for visualisation of missing data,
"DataExplorer", # for visualisation of missing data,
"dplyr",
"magrittr",
"data.table")
# "ivreg") # used to test endogenity redundant
for (i in 1:length(packages)) {
if (!packages[i] %in% rownames(installed.packages())) {
install.packages(packages[i]
, repos = "http://cran.rstudio.com/"
, dependencies = TRUE
)
}
library(packages[i], character.only = TRUE)
}Data
data(CigarettesSW)
stargazer(CigarettesSW,
type = "text", # html, latex
digits = 1)##
## ==============================================================
## Statistic N Mean St. Dev. Min Max
## --------------------------------------------------------------
## cpi 96 1.3 0.2 1.1 1.5
## population 96 5,168,866.0 5,442,345.0 478,447 31,493,524
## packs 96 109.2 25.9 49.3 198.0
## income 96 99,878,736.0 120,541,138.0 6,887,097 771,470,144
## tax 96 42.7 16.1 18.0 99.0
## price 96 143.4 43.9 85.0 240.8
## taxs 96 48.3 19.3 21.3 112.6
## --------------------------------------------------------------CigarettesSW is a panel data set with 96 observations with 9 variables.
- state: factor 48 levels one for each state except Alaska and Hawaii
- year: factor 2 levels 1985 and 1995
- cpi: num consumer price index
- population: num state population
- packs: num number of packs per capita
- income: num state personal income (total)
- tax: num average state,federal, average local excises taxes for fiscal year
- price: num average price during fiscal year(including sales tax)
- taxs: num average excise taxes for fiscal year(including sales tax)
CigarettesSW$ipc <- CigarettesSW$income / CigarettesSW$population
unique(CigarettesSW$cpi) # cpi isn't really a num as cpi doesnt change in year## [1] 1.076 1.524CigarettesSW$cpi <- as.numeric(CigarettesSW$cpi) # creates factor for 2 years
# if there was more years could leave as a num Cleaning
cig_melt <- reshape2::melt(CigarettesSW, id.vars = c("state", "year", "cpi"))
ggplot(data = cig_melt,
mapping = aes(x = value)) +
geom_histogram(bins = 30) +
facet_wrap(. ~ variable, scale = 'free')
# filters to just 1985
cig_1985 <- CigarettesSW %>% filter(year == 1985)
# filter to just 1995
cig_1995 <- CigarettesSW %>% filter(year == 1995)
# create new variables in terms of base year(?) $ cpi 1985
cig_1985_real <- cig_1985 %>%
mutate( # mutate is in 'dplyr' use to create, remove, modifie colms
real_price = price / 1.076,
real_tax = tax / 1.076,
real_taxs = taxs / 1.076,
real_income = income / 1.076,
real_ipc = (income / population) / 1.076
)
# create new variables in terms of base year(?) $ cpi 1995
cig_1995_real <- cig_1995 %>%
mutate(
real_price = price / 1.524,
real_tax = tax / 1.524,
real_taxs = taxs / 1.524,
real_income = income / 1.524,
real_ipc = (income / population) / 1.524
)
# bind the colms back together
CigarettesSW_real <- bind_rows(cig_1985_real, cig_1995_real)cig_melt_real <- reshape2::melt(CigarettesSW_real, id.vars = c("state", "year", "cpi"))
exluded_v <- c("population", "packs", "income")
cig_melt_real %>%
filter(!variable %in% exluded_v) %>%
ggplot(aes(x = value)) +
geom_histogram(bins = 30,) +
facet_wrap(~ variable, scales = "free")
# creates plot order to compare nominal vs real dollars
plot_order <- c(
"price", "real_price",
"tax", "real_tax",
"taxs", "real_taxs",
"ipc", "real_ipc")
cig_melt_real %>%
filter(variable %in% plot_order) %>%
ggplot(aes(x = value)) +
geom_histogram(bins = 30) +
facet_wrap(
~ factor(variable, levels = plot_order),
scales = "free",
ncol = 2
)
1
Is data balanced?
# re-order levels
reorder_size <- function(x) {
factor(x, levels = names(sort(table(x), decreasing = TRUE)))
}
ggplot(data = CigarettesSW_real ,
aes(x = reorder_size(year)
)
) +
geom_bar() +
xlab("Year") + ylab("Frequency") +
theme(axis.text.x = element_text(angle = 45)
)
reorder_size <- function(x) {
factor(x, levels = names(sort(table(x), decreasing = TRUE)))
}
ggplot(data = CigarettesSW_real ,
aes(x = reorder_size(state)
)
) +
geom_bar() +
xlab("state") + ylab("Frequency") +
theme(axis.text.x = element_text(angle = 90)
)
Yes data is balance.
what is the time component and entity component?
Time is year 1985 and 1995. Entity is state.
2 Estimating Equation
The model will look at packs per capita as real price, real IPC, real tax, real taxs(exsise),
\[ \begin{align} \text{Packs}_i \ = &\ \beta_0 \ + \beta_1\text{Price}_i \ + \beta_2\text{IPC}_i \ + \\ &\beta_3\text{Tax}_i \ + \beta_4\text{Taxs}_i \ + \epsilon_i \end{align} \]
\[ \begin{align} \text{Packs}_i \ = &\ \beta_0 \ + \beta_1\text{Real Price}_i \ + \beta_2\text{Real IPC}_i \ + \\ &\beta_3\text{Real Tax}_i \ + \beta_4\text{Real Taxs}_i \ + \epsilon_i \end{align} \]
models
real_ols <- lm(packs ~ real_price + real_ipc + real_tax
+ real_taxs,
data = CigarettesSW_real)
summary(real_ols)##
## Call:
## lm(formula = packs ~ real_price + real_ipc + real_tax + real_taxs,
## data = CigarettesSW_real)
##
## Residuals:
## Min 1Q Median 3Q Max
## -41.880 -11.350 -2.040 9.543 61.165
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 223.7719 16.8479 13.282 <2e-16 ***
## real_price -1.6003 0.2957 -5.413 5e-07 ***
## real_ipc 2.7088 1.1422 2.372 0.0198 *
## real_tax -0.7545 0.9418 -0.801 0.4252
## real_taxs 1.2526 1.0012 1.251 0.2141
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 18.21 on 91 degrees of freedom
## Multiple R-squared: 0.5256, Adjusted R-squared: 0.5048
## F-statistic: 25.21 on 4 and 91 DF, p-value: 4.57e-14vif(real_ols)## real_price real_ipc real_tax real_taxs
## 7.553544 1.712062 19.050219 30.701214There is sever multicollinearity among the predictor variables whether using real dollars or nominal. Second models will be created that will drop tax/taxes as they will be correlated with income per capita \[ \begin{align} \text{Packs}_i \ = &\ \beta_0 \ + \beta_1\text{Real Price}_i \ + \beta_2\text{Real IPC}_i \ + \epsilon_i \end{align} \]
real_ols2 <- lm(packs ~ real_price + real_ipc,
data = CigarettesSW_real)
summary(real_ols2)##
## Call:
## lm(formula = packs ~ real_price + real_ipc, data = CigarettesSW_real)
##
## Residuals:
## Min 1Q Median 3Q Max
## -43.893 -10.057 -1.125 9.860 63.084
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 208.6927 13.4205 15.550 < 2e-16 ***
## real_price -1.2131 0.1341 -9.047 2.12e-14 ***
## real_ipc 2.3055 1.0881 2.119 0.0368 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 18.22 on 93 degrees of freedom
## Multiple R-squared: 0.5142, Adjusted R-squared: 0.5038
## F-statistic: 49.23 on 2 and 93 DF, p-value: 2.622e-15vif(real_ols2)## real_price real_ipc
## 1.550691 1.550691Once taxes have been excluded the model no longer shows multicollinearity. Our R squared has gone down about 1% while adj R squared only went down .1% for the real dollars model.
As the real_ols2 model has the higher R squared and uses real dollars we will move forward with this model. The 208.6927 intercept is the theoretical number of packs sold per capita, if all other variables are held at 0. Since price and income per capital in the real world won’t be 0 this is purely theoretical. The coefficient on price is the expected direction and indicates that a 1 real dollar price increase is associated with a decrease of 1.2131 packs per capita, this was statistically significant at 99.9 confidence interval. The coefficient on real income per capita is the expected sign, it indicates that for a 1 real dollar increase in income per capita it is associated with an increase of 2.3055 packs per capita.This is significant at 95 confidence interval.
The model explains 51.42% of the variance in the packs per capita.
While the real_ols2 models does “correct” the multicollinearity it dropped valuable data from the equation. Dropping the tax variables may cause endogeneity, to test we will perform the Durbin-Wu-Hausman test using instrument variables.
Durbin-Wu-Hausman test
# Run the iv model ivreg runs a 2SLS regression
# replaced ivreg with plm
# pipe to include the iv's
pool_2SLSIV <- plm(packs ~ real_price + real_ipc |
real_tax + real_taxs + real_ipc,
data = CigarettesSW_real,
index = c("state", "year"),
model = "pooling") # pooling turns of FE
summary(pool_2SLSIV, diagnostics = TRUE)## Pooling Model
## Instrumental variable estimation
## (Balestra-Varadharajan-Krishnakumar's transformation)
##
## Call:
## plm(formula = packs ~ real_price + real_ipc | real_tax + real_taxs +
## real_ipc, data = CigarettesSW_real, model = "pooling", index = c("state",
## "year"))
##
## Balanced Panel: n = 48, T = 2, N = 96
##
## Residuals:
## Min. 1st Qu. Median 3rd Qu. Max.
## -45.63615 -10.39839 -0.79849 9.24080 65.46437
##
## Coefficients:
## Estimate Std. Error z-value Pr(>|z|)
## (Intercept) 204.54375 13.75598 14.8694 < 2.2e-16 ***
## real_price -1.11307 0.15087 -7.3779 1.608e-13 ***
## real_ipc 1.82177 1.14035 1.5975 0.1101
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 63586
## Residual Sum of Squares: 31072
## R-Squared: 0.51368
## Adj. R-Squared: 0.50322
## Chisq: 70.9419 on 2 DF, p-value: 3.937e-16summary(pool_2SLSIV)## Pooling Model
## Instrumental variable estimation
## (Balestra-Varadharajan-Krishnakumar's transformation)
##
## Call:
## plm(formula = packs ~ real_price + real_ipc | real_tax + real_taxs +
## real_ipc, data = CigarettesSW_real, model = "pooling", index = c("state",
## "year"))
##
## Balanced Panel: n = 48, T = 2, N = 96
##
## Residuals:
## Min. 1st Qu. Median 3rd Qu. Max.
## -45.63615 -10.39839 -0.79849 9.24080 65.46437
##
## Coefficients:
## Estimate Std. Error z-value Pr(>|z|)
## (Intercept) 204.54375 13.75598 14.8694 < 2.2e-16 ***
## real_price -1.11307 0.15087 -7.3779 1.608e-13 ***
## real_ipc 1.82177 1.14035 1.5975 0.1101
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 63586
## Residual Sum of Squares: 31072
## R-Squared: 0.51368
## Adj. R-Squared: 0.50322
## Chisq: 70.9419 on 2 DF, p-value: 3.937e-16stargazer(real_ols2, pool_2SLSIV,
type = "text",
omit = c("state", "year"),
#column.labels = c("real ols2", "real iv"),
notes = "dropping state and year coef")##
## =====================================================
## Dependent variable:
## ---------------------------------
## packs
## OLS panel
## linear
## (1) (2)
## -----------------------------------------------------
## real_price -1.213*** -1.113***
## (0.134) (0.151)
##
## real_ipc 2.305** 1.822
## (1.088) (1.140)
##
## Constant 208.693*** 204.544***
## (13.421) (13.756)
##
## -----------------------------------------------------
## Observations 96 96
## R2 0.514 0.514
## Adjusted R2 0.504 0.503
## Residual Std. Error 18.224 (df = 93)
## F Statistic 49.226*** (df = 2; 93) 70.942***
## =====================================================
## Note: *p<0.1; **p<0.05; ***p<0.01
## dropping state and year coefWith a p-value of 3.559e-12 we reject the null hypothesis and will use a 2SLS model which will take price as a function of tax + taxes.
2 Stage Equation
\[ \begin{align} \text{First stage:} \ \text{Real Price}_{st} \ = &\ \gamma_s \ + \delta_1\text{Real Tax}_{st} \ + \delta_2\text{Real Taxs}_{st} \ + \delta_3\text{Real IPC}_{st} \ + \nu_{st} \\ \text{Second stage:} \ \text{Packs}_{st} \ = &\ \alpha_s \ + \beta_1\widehat{\text{Price}_{st}} \ + \beta_2\text{Real IPC}_{st} \ + \epsilon_{st} \end{align} \]
Where \(\gamma_s\) is the state fixed effects on the first stage and \(\nu_{st}\) is the first stage error term.
Where s represents states, t represent time in year. \(\alpha_s\) holds the constant effects of the states that do not change over time, while \(\epsilon_{st}\) captures the non-constant effects over time. Where \(\widehat{\text{Price}_{st}}\) is Real Price
3 Fixed Effects
\[ \begin{align} \text{Packs}_{st} \ = &\ \alpha_s \ + \delta_t + \beta_1\text{Real Price}_{st} \ + \beta_2\text{Real IPC}_{st} \ + \epsilon_{st} \end{align} \]
Where s represents states, t represent time in year. \(\alpha_s\) holds the constant effects of the states that do not change over time, and \(\delta_t\) represents the constant effects over time, while \(\epsilon_{st}\) captures the non-constant effects over time.
Method 1 : Using Two Way Fixed Effects
Using the panel command we can expand the model by adding in fixed effects for year and state.
# pipe to include the iv's
TWFEIV <- plm(packs ~ real_price + real_ipc |
real_tax + real_taxs + real_ipc,
data = CigarettesSW_real,
index = c("state", "year"),
model = "within",
effect = "twoways")
# within puts the FE into model
# twoways says both year and state fe
stargazer(TWFEIV,
type = "text",
omit = c("state", "year"))##
## ========================================
## Dependent variable:
## ---------------------------
## packs
## ----------------------------------------
## real_price -0.779***
## (0.114)
##
## real_ipc -0.798
## (1.840)
##
## ----------------------------------------
## Observations 96
## R2 0.556
## Adjusted R2 0.063
## F Statistic 48.874***
## ========================================
## Note: *p<0.1; **p<0.05; ***p<0.01Note that the adj \(r^2\) is 0.063 this can be ignored because of the way stargazer calculates \(r^2\) it takes each year, and state as a new variable artificially lowering adj \(r^2\). We will remove from our comparison.
stargazer(pool_2SLSIV, TWFEIV,
type = "text",
omit = c("state", "year"),
column.labels = c("pooled 2SLSIV", "TWFEIV"),
notes = "dropping state and year coef",
omit.stat = "adj.rsq")##
## ==========================================
## Dependent variable:
## -----------------------------
## packs
## pooled 2SLSIV TWFEIV
## (1) (2)
## ------------------------------------------
## real_price -1.113*** -0.779***
## (0.151) (0.114)
##
## real_ipc 1.822 -0.798
## (1.140) (1.840)
##
## Constant 204.544***
## (13.756)
##
## ------------------------------------------
## Observations 96 96
## R2 0.514 0.556
## F Statistic 70.942*** 48.874***
## ==========================================
## Note: *p<0.1; **p<0.05; ***p<0.01
## dropping state and year coefThis shows that by accounting for the FE of state and time strips away some of the error for price making the model more precise.
The sign on real_ipc coefficient flipped to negative. The reason for the flip is because the FE are now controlling more accurately for increases in income. Cigarettes are highly addictive, in both the pooled 2SLSIV and TWFEIV is not statistically significant meaning they are indistinguishable from 0, this along with the sign flip shows that the cigarettes are demand inelastic cigarette consumers will still continue to buy cigarettes at roughly the same rate whether a persons real income increase, due time or across state lines.
2 2SLS IV model with dummies
\[ \begin{align} \text{First stage:} \ \text{Real Price}_{i} \ = &\ \delta_0 \ + \delta_1\text{Real Tax}_{i} \ + \delta_2\text{Real Taxs}_{i} \ + \delta_3\text{Real IPC}_{i} \ + \\ &\ \delta_4\text{state}_{i} \ + \delta_5\text{year}_{i} + \nu_{i} \\ \text{Second stage:} \ \text{Packs}_{i} \ = &\ \beta_0 \ + \beta_1\widehat{\text{Price}_{i}} \ + \beta_2\text{Real IPC}_{i} \ + \beta_3\text{State}_{i} + \beta_4\text{year}_{i} + \epsilon_{i} \end{align} \] Since my base model is not an ols, and includes IV, we will have to make use some weird code to model this.
# Instead of let plm create dummies we manually put them into the model
# this essentially acts as if we did a normal OLS w/ FE
FE_2SLSIV <- plm(packs ~ real_price + real_ipc + factor(state) + factor(year) |
real_tax + real_taxs + real_ipc + factor(state) + factor(year),
data = CigarettesSW_real,
index = c("state", "year"),
model = "pooling") # pooling turns of FE
stargazer(FE_2SLSIV, TWFEIV,
type = "text",
omit = c("state", "year"),
column.labels = c("FE 2SLSIV", "TWFEIV"),
notes = "dropping state and year coef",
omit.stat = "adj.rsq")##
## ==========================================
## Dependent variable:
## -----------------------------
## packs
## FE 2SLSIV TWFEIV
## (1) (2)
## ------------------------------------------
## real_price -0.779*** -0.779***
## (0.114) (0.114)
##
## real_ipc -0.798 -0.798
## (1.840) (1.840)
##
## Constant 198.550***
## (21.591)
##
## ------------------------------------------
## Observations 96 96
## R2 0.976 0.556
## F Statistic 1,838.281*** 48.874***
## ==========================================
## Note: *p<0.1; **p<0.05; ***p<0.01
## dropping state and year coefUsing the LSDV method applied to my 2SLSIV by adding in the dummy variables we were able to show that the the coefficients remain the same.
3 Demeaned Regression
\[ \begin{align} \text{First stage:} \ \widetilde{\text{Real Price}}_{st} \ = &\ \delta_1\widetilde{\text{Real Tax}}_{st} \ + \ \delta_2\widetilde{\text{Real Taxs}}_{st} \ + \ \delta_3\widetilde{\text{Real IPC}}_{st} \ + \ \nu_{st} \\ \text{Second stage:} \ \widetilde{\text{Packs}}_{st} \ = &\ \beta_1\widehat{\widetilde{\text{Price}}}_{st} \ + \ \beta_2\widetilde{\text{Real IPC}}_{st} \ + \ \epsilon_{st} \end{align} \]
# Properly create the two-way demeaned data frame (including packs!)
cig_demeaned <- data.frame(
packs = CigarettesSW_real$packs -
ave(CigarettesSW_real$packs, CigarettesSW_real$state) -
ave(CigarettesSW_real$packs, CigarettesSW_real$year) +
mean(CigarettesSW_real$packs),
real_price = CigarettesSW_real$real_price -
ave(CigarettesSW_real$real_price, CigarettesSW_real$state) -
ave(CigarettesSW_real$real_price, CigarettesSW_real$year) +
mean(CigarettesSW_real$real_price),
real_ipc = CigarettesSW_real$real_ipc -
ave(CigarettesSW_real$real_ipc, CigarettesSW_real$state) -
ave(CigarettesSW_real$real_ipc, CigarettesSW_real$year) +
mean(CigarettesSW_real$real_ipc),
real_tax = CigarettesSW_real$real_tax -
ave(CigarettesSW_real$real_tax, CigarettesSW_real$state) -
ave(CigarettesSW_real$real_tax, CigarettesSW_real$year) +
mean(CigarettesSW_real$real_tax),
real_taxs = CigarettesSW_real$real_taxs -
ave(CigarettesSW_real$real_taxs, CigarettesSW_real$state) -
ave(CigarettesSW_real$real_taxs, CigarettesSW_real$year) +
mean(CigarettesSW_real$real_taxs),
state = CigarettesSW_real$state,
year = CigarettesSW_real$year
)
# we will again run our model through a modified plm with manual dummy variable
demeaned_model <-
plm(packs ~ real_price + real_ipc + factor(state)
+ factor(year) | real_tax + real_taxs + real_ipc
+ factor(state) + factor(year),
data = cig_demeaned,
index = c("state", "year"),
model = "pooling")
# View results
stargazer(demeaned_model,
type = "text",
omit = c("state", "year"))##
## ========================================
## Dependent variable:
## ---------------------------
## packs
## ----------------------------------------
## real_price -0.779***
## (0.114)
##
## real_ipc -0.798
## (1.840)
##
## Constant -0.000
## (4.143)
##
## ----------------------------------------
## Observations 96
## R2 0.556
## Adjusted R2 0.063
## F Statistic 48.874
## ========================================
## Note: *p<0.1; **p<0.05; ***p<0.01The demeaned section code I had to used ai to build out the variables because of having two way fixed effects, I was unsure how to structure the math to get the correct means. The set up is to take the average of variable with respect to state, subtract the average of the variable with respect to year add the mean of the variable. this is then subtracted from the variable itself. What this does it takes the states fixed effects, the years fixed effects and strips it out of the variable. This then causes the same Fixed Effects calculation which results in the same coefficients given by the previous two models.
stargazer( TWFEIV, FE_2SLSIV, demeaned_model,
type = "text",
omit = c("state", "year"),
column.labels = c( "TWFEIV","FE2SLSIV", "Demeaned"),
notes = "dropping state and year coef",
omit.stat = "adj.rsq")##
## =============================================
## Dependent variable:
## --------------------------------
## packs
## TWFEIV FE2SLSIV Demeaned
## (1) (2) (3)
## ---------------------------------------------
## real_price -0.779*** -0.779*** -0.779***
## (0.114) (0.114) (0.114)
##
## real_ipc -0.798 -0.798 -0.798
## (1.840) (1.840) (1.840)
##
## Constant 198.550*** -0.000
## (21.591) (4.143)
##
## ---------------------------------------------
## Observations 96 96 96
## R2 0.556 0.976 0.556
## F Statistic 48.874*** 1,838.281*** 48.874
## =============================================
## Note: *p<0.1; **p<0.05; ***p<0.01
## dropping state and year coefThe Coefficients do not change because each of the three methods is essentially doing the same thing. Each model is using a different way to strip the variance from state and time out of the model. They all controll for these fixed effects in their own way.
It is common to use both time and entity when using palen data. Panel data is of a certain group, in my case states, and over time, years in my case, since panel data includes both it often makes sense to control for both when creating models.
The three models each specify the fixed effects in their own manner and the results show that the coefficients do not change.