Cigarette Consumption Data

This is the last lab, which is largely a review of lab 6 but this time we will be using Chris’s code from topic 9. See below the codebook of the data.

Variable Description1
cpi consumer price index
pop state population
packpc number of packs per capita
income95 state personal income (total)
income95pc state personal income (per capita)
tax95 average state, federal, and average local excise taxes for fiscal year
taxs95 average excise taxes for fiscal year, including sales taxes
avgprs95 average price during fiscal year, including sales taxes (cents)
pretax95 pretax95 = avgprs95 - taxs95
# Load cigarette consumption data (Jonathan Gruber, MIT)
# Variables (see codebook):
# state year    cpi pop packpc  income  tax avgprs  taxs
df <- read.csv("https://faculty.washington.edu/cadolph/panUW/cigarette.csv")

# Quick inflation adjustment to 1995 dollars
inflAdjust <- function(x,cpi,year,target) {
  unique(cpi[year==target])*x/cpi
}


df <-
  df %>% 
  mutate(
    # adjust series for inflation
    income95 = inflAdjust(income, cpi, year, 1995),
    tax95 = inflAdjust(tax, cpi, year, 1995),
    avgprs95 = inflAdjust(avgprs, cpi, year, 1995),
    taxs95 = inflAdjust(taxs, cpi, year, 1995),
    # create a per capita income (in k)
    income95pc = income95/pop,
    # Create pretax price, 1995 dollars
    pretax95 = avgprs95 - taxs95
  )



## Construct the weighted avg of consumption over states
(stateList <- unique(as.character(df$state)))
##  [1] "AL" "AR" "AZ" "CA" "CO" "CT" "DE" "FL" "GA" "IA" "ID" "IL" "IN" "KS" "KY"
## [16] "LA" "MA" "MD" "ME" "MI" "MN" "MO" "MS" "MT" "NC" "ND" "NE" "NH" "NJ" "NM"
## [31] "NV" "NY" "OH" "OK" "OR" "PA" "RI" "SC" "SD" "TN" "TX" "UT" "VA" "VT" "WA"
## [46] "WI" "WV" "WY"
(nstate <- length(stateList))
## [1] 48
(yearList <- unique(df$year))
##  [1] 1985 1986 1987 1988 1989 1990 1991 1992 1993 1994 1995
(nyear <- length(yearList))
## [1] 11
packpcMat <- taxs95Mat <- popMat <- NULL


for (istate in 1:nstate) {
  # for the istate
  curstate <- stateList[istate]
  # subset their packpc and bind into vector
  packpcMat <- cbind(packpcMat, df$packpc[as.character(df$state)==curstate])
  # subset their taxs95 and bind into vector
  taxs95Mat <- cbind(taxs95Mat, df$taxs95[as.character(df$state)==curstate])
  # subset their pop and bind into vector
  popMat <- cbind(popMat, df$pop[as.character(df$state)==curstate])
}

## Take the weighted average of packpc and tax 
(packpcWAvg <- apply(packpcMat*popMat, 1, sum)/apply(popMat, 1, sum))
##  [1] 119.82600 117.01874 113.95779 110.77500 106.27144 100.56928  96.24902
##  [8]  93.44560  90.88317  88.57571  88.69922
(taxs95WAvg <- apply(taxs95Mat*popMat, 1, sum)/apply(popMat, 1, sum))
##  [1] 50.00587 50.46155 50.35367 50.90960 51.94096 54.06268 57.85670 61.18315
##  [9] 62.42414 66.23973 67.35962
# Make a list of complete state names using built in objects
stateFull <- sapply(stateList, function(x) {state.name[grep(x, state.abb)]})

# First, create a plm data frame (special data frame that "knows" the
# unit variable and time variable
pdata <- pdata.frame(df, index=c("state", "year"))

Estimate Arellano-Bond GMM for fixed effects with lagged DV

Remember, pgmm needs formulas in a specific format: 1. in the first part of the RHS, include lags of DV and covariates, as shown 2. in the second part, include the panel data instruments (99 here means use up to the 99th lag of the difference as an instrument) 3. in an optional (not shown) third part of the RHS, include any other instruments

Sargan test has a null of the instruments as a group being exogenous- higher p is good. The residuals of the differences equations should exhibit AR(1) but not AR(2) behavior higher p for AR(2) test is good

# note that pgmm formulas construct lag() properly for panel data,
# though lag() usually doesn't
pgmmformula.3a <- log(packpc) ~ lag(log(packpc), 1) + log(income95pc) + log(avgprs95) | lag(log(packpc), 2:99)


# Model 3: Elasticity specification

# Try difference GMM with two way effects
pgmm.res3g <- pgmm(pgmmformula.3a,
                   data = pdata,
                   effect = "twoways",   # should consider two-way for small T
                   transformation = "d")    # should do ld if T=3

summary(pgmm.res3g)
## Twoways effects One-step model Difference GMM 
## 
## Call:
## pgmm(formula = pgmmformula.3a, data = pdata, effect = "twoways", 
##     transformation = "d")
## 
## Balanced Panel: n = 48, T = 11, N = 528
## 
## Number of Observations Used: 432
## Residuals:
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## -0.181013 -0.019492 -0.001678  0.000000  0.017812  0.204343 
## 
## Coefficients:
##                     Estimate Std. Error z-value  Pr(>|z|)    
## lag(log(packpc), 1)  0.33315    0.14500  2.2976   0.02158 *  
## log(income95pc)      0.18131    0.18166  0.9981   0.31825    
## log(avgprs95)       -0.62271    0.11127 -5.5965 2.188e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Sargan test: chisq(44) = 43.30027 (p-value = 0.5015)
## Autocorrelation test (1): normal = -3.620058 (p-value = 0.00029454)
## Autocorrelation test (2): normal = 0.5219654 (p-value = 0.60169)
## Wald test for coefficients: chisq(3) = 45.27325 (p-value = 8.0945e-10)
## Wald test for time dummies: chisq(9) = 71.0946 (p-value = 9.2856e-12)
# Good Sargan test, Good AR(2) test, Wald supports 2-way

Simulate conditional forecasts from model 3g

# Forecast for 3 years from 1996 to 1998
periods.out <- 3
sims <- 1000

# How big a change in price to simulate?
# How about "double" the average tax in the most recent year?
summary(pdata$taxs95[pdata$year==1995])
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   34.44   48.75   59.84   61.87   74.78  112.63
# The average (and median) tax is about 60 cents/pack
sd(pdata$taxs95[pdata$year==1995])
## [1] 18.47741
# A 60 cent increase would also be about 3 sd's,
# and raise the tax to a bit more than the max observed

# Other possibilities:
# (2) A 10 cent increase
# (3) Raise every sate to the max observed for any state in 1995 (112.60 cents)


# Construct matrix with year dummies
yearfe <- makeFEdummies(pdata$year)
yearfe <- yearfe[,3:ncol(yearfe)]  # Why drop first 2 col's?

# create a formula with years-dummies
yearlist <- unique(pdata$year)
yearlist <- yearlist[3:length(yearlist)]
colnames(yearfe) <- paste0("y",yearlist)

# Construct formulas -- with year dummies (1g)
formula <- "packpc ~ income95pc + avgprs95 -1"

# Bind together the pdata with the matrix of dummies
datayearfe <- cbind(pdata,yearfe)

yearfenames <- NULL # create an empty vector

# loop the formula to incorporate all the years
for (i in 1:ncol(yearfe)) {
  formula <- paste0(formula,"+ y",yearlist[i]," ")
  yearfenames <- c(yearfenames,paste0("y",yearlist[i]))
}

names(datayearfe) <- c(names(df),yearfenames)

formula.1g <- as.formula(formula)

Forecast: average state raised taxes 60 cents

Make matrix of hypothetical x’s:

  • Assume an average state raised taxes 60 cents starting 1996.

Issues: we need to somehow include the state and year FEs. Let’s set the state to be an “average” state in 1995, and year to be like the last year (1995)

pgmm uses covariates in differenced form so we want most of them to be 0 (no change) exceptions: 1 - changes in covariates of interest 2 - time dummies aren’t differenced

# Forecast: Model 3g, +60

simparam.3g <- mvrnorm(sims, coefficients(pgmm.res3g), vcovHC(pgmm.res3g))
simphis.3g <- simparam.3g[,1]
simbetas.3g <- simparam.3g[,2:ncol(simparam.3g)]

# Make matrix of hypothetical x's: covariates
# Still use the 1g formula (no logs) -- we will handle logging manually
#  to get the differences of logs right
xhyp.3g <- cfMake(formula = formula.1g, 
                  data = datayearfe, 
                  nscen = periods.out)

# pgmm uses covariates in differenced form
# so we want most of them to be 0 (no change)
changeTax <- c(60,0,0)

# Need log version of differenced key covariate (doubling tax in avg state)
(meanPrice <- mean(pdata$avgprs, na.rm=TRUE))
## [1] 150.2528
(meanIncome <- mean(pdata$income95pc, na.rm=TRUE))
## [1] 21.25106
xhyp.3g <- cfChange(xscen = xhyp.3g, 
                    covname = "avgprs95",
                    x=log(meanPrice+changeTax[1]) - log(meanPrice),
                    xpre=log(meanPrice) - log(meanPrice),
                    scen=1)

xhyp.3g <- cfChange(xscen = xhyp.3g,  
                    covname = "avgprs95",
                    x=log(meanPrice+changeTax[2]) - log(meanPrice),
                    xpre=log(meanPrice) - log(meanPrice),
                    scen=2)

xhyp.3g <- cfChange(xscen = xhyp.3g,  
                    covname = "avgprs95",
                    x=log(meanPrice+changeTax[3]) - log(meanPrice),
                    xpre=log(meanPrice) - log(meanPrice),
                    scen=3)

xhyp.3g <- cfChange(xscen = xhyp.3g,  
                    covname = "income95pc",
                    x=log(meanIncome) - log(meanIncome),
                    xpre=log(meanIncome) - log(meanIncome),
                    scen=1:3)

for (iyear in 1985:1991){
  xhyp.3g <- cfChange(xscen = xhyp.3g, 
                      covname = paste0("y",iyear), 
                      x=0, 
                      xpre=0, 
                      scen=1:3)
}

xhyp.3g <- cfChange(xscen = xhyp.3g, 
                    covname = "y1992", x=-1, xpre=-1, scen=1)
xhyp.3g <- cfChange(xscen = xhyp.3g,  
                    covname = "y1992", x=0, xpre=0, scen=2)
xhyp.3g <- cfChange(xscen = xhyp.3g,  
                    covname = "y1992", x=0, xpre=0, scen=3)

xhyp.3g <- cfChange(xscen = xhyp.3g,  
                    covname = "y1993", x=1, xpre=1, scen=1)
xhyp.3g <- cfChange(xscen = xhyp.3g,  
                    covname = "y1993", x=-1, xpre=-1, scen=2)
xhyp.3g <- cfChange(xscen = xhyp.3g,  
                    covname = "y1993", x=0, xpre=0, scen=3)

xhyp.3g <- cfChange(xscen = xhyp.3g,  
                    covname = "y1994", x=0, xpre=0, scen=1)
xhyp.3g <- cfChange(xscen = xhyp.3g,  
                    covname = "y1994", x=1, xpre=1, scen=2)
xhyp.3g <- cfChange(xscen = xhyp.3g,  
                    covname = "y1994", x=-1, xpre=-1, scen=3)

xhyp.3g <- cfChange(xscen = xhyp.3g,  
                    covname = "y1995", x=0, xpre=0, scen=1)
xhyp.3g <- cfChange(xscen = xhyp.3g,  
                    covname = "y1995", x=0, xpre=0, scen=2)
xhyp.3g <- cfChange(xscen = xhyp.3g,  
                    covname = "y1995", x=1, xpre=1, scen=3)

# We can "ignore" the state fixed effects for now and add them later
# because model is total linear

# Create baseline scenario
xbase.3g <- xhyp.3g
xbase.3g$x <- xbase.3g$xpre

# We need a lag of the price per pack
lagY.3g <- NULL # Hypothetical previous change in Y for simulation

for (iunit in 1:length(pgmm.res3g$model)) {
  lagY.3g <- c(lagY.3g, as.data.frame(pgmm.res3g$model[[iunit]])["1995",1])
}

lagY.3g <- mean(lagY.3g, na.rm=TRUE)


initialY <- mean(pdata$packpc[pdata$year==1993], na.rm=TRUE)
# Hypothetical initial level of Y for simulation

# Simulate expected values of Y (on original level scale)
# out to periods.out given hypothetical future values of X,
# initial lags of the change in Y, and an initial level of Y
sim.ev3g <- ldvsimev(xhyp.3g,               # The matrix of hypothetical x's
                     simbetas.3g,           # The matrix of simulated betas
                     ci=0.95,            # Desired confidence interval
                     constant=NA,        # NA indicates no constant!
                     phi=simphis.3g,            # estimated AR parameters; length must match lagY 
                     lagY=lagY.3g,          # lags of y, most recent last
                     transform="difflog",   # "log" to undo log transformation,
                     # "diff" to under first differencing
                     # "difflog" to do both
                     initialY=initialY   # for differenced models, the lag of the level of y
)

# Simulate expected values of Y given no change in covariates
sim.base3g <- ldvsimev(xbase.3g,               # The matrix of hypothetical x's
                       simbetas.3g,           # The matrix of simulated betas
                       ci=0.95,            # Desired confidence interval
                       constant=NA,        # NA indicates no constant!
                       phi=simphis.3g,            # estimated AR parameters; length must match lagY 
                       lagY=lagY.3g,          # lags of y, most recent last
                       transform="difflog",   # "log" to undo log transformation,
                       # "diff" to under first differencing
                       # "difflog" to do both
                       initialY=initialY   # for differenced models, the lag of the level of y
) 

# Simulate first differences in y (as cumulated changes)
# out to periods.out given hypothetical future values of x, xpre,
# and initial lags of the change in y
sim.fd3g <- ldvsimfd(xhyp.3g,            # The matrix of hypothetical x's
                     simbetas.3g,        # The matrix of simulated betas
                     ci=0.95,            # Desired confidence interval
                     constant=NA,        # Column containing the constant
                     # set to NA for no constant
                     phi=simphis.3g,     # estimated AR parameters; length must match lagY 
                     lagY=lagY.3g,       # lags of y, most recent last
                     transform="difflog",# Model is differenced logs
                     initialY=initialY   # Required in this case (fd of differenced log Y)
)

# Simulate relative risks in y
# out to periods.out given hypothetical future values of x, xpre,
# and initial lags of the change in y
sim.rr3g <- ldvsimrr(xhyp.3g,            # The matrix of hypothetical x's
                     simbetas.3g,        # The matrix of simulated betas
                     ci=0.95,            # Desired confidence interval
                     constant=NA,        # Column containing the constant
                     # set to NA for no constant
                     phi=simphis.3g,     # estimated AR parameters; length must match lagY 
                     lagY=lagY.3g,       # lags of y, most recent last
                     transform="difflog",# Model is differenced logs
                     initialY=initialY   # Required for differenced Y in ldvsimrr
)

Forecast: avg state, by state, and weighted avg of all states

In this part we will repeat the same program as before, but we will loop the code such that we can estimate the impact in each state, and also take the weighted average of all the states.

Quantities of interest: - packs FD - packs RR/%change - revenue pc FD - revenue mils FD - revenue %GSP FD

The loop below will compute all these quantities for each state. Note that we will be extracting all the quantities using Chris’s custom functions collectUnitSims and aggregateUnitSims. These functions are included in the the TS_code.R script, in the source folder.

# Visuals: line plots.  1x2 and 1x3 for two diff scenarios (x2)

# List for results
simEV.3g.scen1 <- simBASE.3g.scen1 <- simFD.3g.scen1 <- simRR.3g.scen1 <- list()

# All states modeled
allstates <- names(pgmm.res3g$model)

# Storage vector for weights
popweight <- rep(NA, length(allstates))

# Loop over states in model
for (i in 1:length(pgmm.res3g$model)) {
  
  curstate <- allstates[i]
  
  # Grab population weight for later
  popweight[i] <- pdata$pop[(pdata$year==1995)&(pdata$state==curstate)]
  
  # Make matrix of hypothetical x's: covariates
  # Still use the 1g formula (no logs) -- we will handle logging manually
  #  to get the differences of logs right
  xhyp.3g <- cfMake(formula.1g, datayearfe, periods.out)
  
  # pgmm uses covariates in differenced form
  # so we want most of them to be 0 (no change)
  # exceptions:
  # (1) changes in covariates of interest
  # (2) time dummies aren't differenced
  #xhyp.3g$x <- xhyp.3g$xpre <- 0*xhyp.3g$x
  
  # Need log version of differenced key covariate (doubling tax in avg state)
  curPrice92 <- pdata$avgprs95[(pdata$state==curstate)&(pdata$year==1992)]
  curPrice93 <- pdata$avgprs95[(pdata$state==curstate)&(pdata$year==1993)]
  curPrice94 <- pdata$avgprs95[(pdata$state==curstate)&(pdata$year==1994)]
  curPrice95 <- pdata$avgprs95[(pdata$state==curstate)&(pdata$year==1995)]
  curIncome92 <-  pdata$income95pc[(pdata$state==curstate)&(pdata$year==1992)]
  curIncome93 <-  pdata$income95pc[(pdata$state==curstate)&(pdata$year==1993)]
  curIncome94 <-  pdata$income95pc[(pdata$state==curstate)&(pdata$year==1994)]
  curIncome95 <-  pdata$income95pc[(pdata$state==curstate)&(pdata$year==1995)]
  
  changeTax <- c(60, 0, 0)
  
  xhyp.3g <- cfChange(xhyp.3g, "avgprs95",
                      x=log(curPrice93+changeTax[1]) - log(curPrice92),
                      xpre=log(curPrice93) - log(curPrice92),
                      scen=1)
  
  xhyp.3g <- cfChange(xhyp.3g, "avgprs95",
                      x=log(curPrice94+changeTax[2]) - log(curPrice93),
                      xpre=log(curPrice94) - log(curPrice93),
                      scen=2)
  
  xhyp.3g <- cfChange(xhyp.3g, "avgprs95",
                      x=log(curPrice95+changeTax[3]) - log(curPrice94),
                      xpre=log(curPrice95) - log(curPrice94),
                      scen=3)
  
  xhyp.3g <- cfChange(xhyp.3g, "income95pc",
                      x=log(curIncome93) - log(curIncome92),
                      xpre=log(curIncome93) - log(curIncome92),
                      scen=1)
  
  xhyp.3g <- cfChange(xhyp.3g, "income95pc",
                      x=log(curIncome94) - log(curIncome93),
                      xpre=log(curIncome94) - log(curIncome93),
                      scen=2)
  
  xhyp.3g <- cfChange(xhyp.3g, "income95pc",
                      x=log(curIncome95) - log(curIncome94),
                      xpre=log(curIncome95) - log(curIncome94),
                      scen=3)
  
  for (iyear in 1985:1991){
    xhyp.3g <- cfChange(xhyp.3g, paste0("y",iyear), x=0, xpre=0, scen=1:3)
  }
  
  xhyp.3g <- cfChange(xhyp.3g, "y1992", x=-1, xpre=-1, scen=1)
  xhyp.3g <- cfChange(xhyp.3g, "y1992", x=0, xpre=0, scen=2)
  xhyp.3g <- cfChange(xhyp.3g, "y1992", x=0, xpre=0, scen=3)
  
  xhyp.3g <- cfChange(xhyp.3g, "y1993", x=1, xpre=1, scen=1)
  xhyp.3g <- cfChange(xhyp.3g, "y1993", x=-1, xpre=-1, scen=2)
  xhyp.3g <- cfChange(xhyp.3g, "y1993", x=0, xpre=0, scen=3)
  
  xhyp.3g <- cfChange(xhyp.3g, "y1994", x=0, xpre=0, scen=1)
  xhyp.3g <- cfChange(xhyp.3g, "y1994", x=1, xpre=1, scen=2)
  xhyp.3g <- cfChange(xhyp.3g, "y1994", x=-1, xpre=-1, scen=3)
  
  xhyp.3g <- cfChange(xhyp.3g, "y1995", x=0, xpre=0, scen=1)
  xhyp.3g <- cfChange(xhyp.3g, "y1995", x=0, xpre=0, scen=2)
  xhyp.3g <- cfChange(xhyp.3g, "y1995", x=1, xpre=1, scen=3)
  
  # Baseline scenario
  xbase.3g <- xhyp.3g
  xbase.3g$x <- xbase.3g$xpre
  
  # We need a lag of the price per pack
  lagY.3g <- as.data.frame(pgmm.res3g$model[[i]])["1992",1]
  
  # This restores the fixed effect by state
  initialY <- pdata$packpc[(pdata$state==curstate)&(pdata$year==1992)]
  
  
  # Simulate expected values of Y (on original level scale)
  # out to periods.out given hypothetical future values of X,
  # initial lags of the change in Y, and an initial level of Y
  simState.ev3g <- ldvsimev(xhyp.3g,               # The matrix of hypothetical x's
                            simbetas.3g,           # The matrix of simulated betas
                            ci=0.95,            # Desired confidence interval
                            constant=NA,        # NA indicates no constant!
                            phi=simphis.3g,            # estimated AR parameters; length must match lagY 
                            lagY=lagY.3g,          # lags of y, most recent last
                            transform="difflog",   # "log" to undo log transformation,
                            # "diff" to under first differencing
                            # "difflog" to do both
                            initialY=initialY,   # for differenced models, the lag of the level of y
                            simulates=TRUE
  )
  
  # Simulate expected values of Y given no change in covariates
  simState.base3g <- ldvsimev(xbase.3g,               # The matrix of hypothetical x's
                              simbetas.3g,           # The matrix of simulated betas
                              ci=0.95,            # Desired confidence interval
                              constant=NA,        # NA indicates no constant!
                              phi=simphis.3g,            # estimated AR parameters; length must match lagY 
                              lagY=lagY.3g,          # lags of y, most recent last
                              transform="difflog",   # "log" to undo log transformation,
                              # "diff" to under first differencing
                              # "difflog" to do both
                              initialY=initialY,   # for differenced models, the lag of the level of y
                              simulates=TRUE
  ) 
  
  # Simulate first differences in y (as cumulated changes)
  # out to periods.out given hypothetical future values of x, xpre,
  # and initial lags of the change in y
  simState.fd3g <- ldvsimfd(xhyp.3g,            # The matrix of hypothetical x's
                            simbetas.3g,        # The matrix of simulated betas
                            ci=0.95,            # Desired confidence interval
                            constant=NA,        # Column containing the constant
                            # set to NA for no constant
                            phi=simphis.3g,     # estimated AR parameters; length must match lagY 
                            lagY=lagY.3g,       # lags of y, most recent last
                            transform="difflog",# Model is differenced logs
                            initialY=initialY,  # Required in this case (fd of differenced log Y)
                            simulates=TRUE
  )
  
  # Simulate relative risks in y
  # out to periods.out given hypothetical future values of x, xpre,
  # and initial lags of the change in y
  simState.rr3g <- ldvsimrr(xhyp.3g,            # The matrix of hypothetical x's
                            simbetas.3g,        # The matrix of simulated betas
                            ci=0.95,            # Desired confidence interval
                            constant=NA,        # Column containing the constant
                            # set to NA for no constant
                            phi=simphis.3g,     # estimated AR parameters; length must match lagY 
                            lagY=lagY.3g,       # lags of y, most recent last
                            transform="difflog",# Model is differenced logs
                            initialY=initialY,  # Required for differenced Y in ldvsimrr
                            simulates=TRUE
  )
  
  
  # Collect results (only works for 1 CI)
  simEV.3g.scen1 <- collectUnitSims(simEV.3g.scen1, simState.ev3g)
  simBASE.3g.scen1 <- collectUnitSims(simBASE.3g.scen1, simState.base3g)
  simFD.3g.scen1 <- collectUnitSims(simFD.3g.scen1, simState.fd3g)
  simRR.3g.scen1 <- collectUnitSims(simRR.3g.scen1, simState.rr3g)
}

# Compute weighted average across states
simEV.3g.scen1 <- aggregateUnitSims(simEV.3g.scen1, weighted.mean, w=popweight)
simBASE.3g.scen1 <- aggregateUnitSims(simBASE.3g.scen1, weighted.mean, w=popweight)
simFD.3g.scen1 <- aggregateUnitSims(simFD.3g.scen1, weighted.mean, w=popweight)
simRR.3g.scen1 <- aggregateUnitSims(simRR.3g.scen1, weighted.mean, w=popweight)

Plot the FD and RR results using a custom tile-graphics function

Finally, we will plot all the quantities of interest using Chris’s tile package.

# Nice colors
col <- brewer.pal(4,"Dark2")[3:4]

# Other settings common to all plots
avg <- "weighted.mean"
periods <- 1993:1995
limits <- c(1991.9, 1995.25, -43, 2)
xat <- c(1992, 1993, 1994, 1995)
plottitles <- c("Change, Packs pc after +$.60 tax",
                "Percent Change, Packs pc after +$.60 tax")
xtitles <- c("Forecast Year", "Forecast Year")
ytitles <- c("","")

# Plot the results by state
cigUnitLineplots(FD=simFD.3g.scen1, RR=simRR.3g.scen1,
                 periods=periods, units=stateList,
                 showUnits=TRUE, showMean=FALSE,
                 avg=avg, avgLab="US",
                 unitLabX=1995.125, labSet="tb", labK=1,
                 initialZero=TRUE, CI=TRUE,
                 limitsFD=limits, limitsRR=limits, xat=xat,
                 col=col, dots=TRUE, RRasPD=TRUE,
                 plottitles=plottitles,
                 xtitles=xtitles, ytitles=ytitles)

# Plot the results by state with weighted mean
cigUnitLineplots(FD=simFD.3g.scen1, RR=simRR.3g.scen1,
                 periods=periods, units=stateList,
                 showUnits=TRUE, showMean=TRUE,
                 avg=avg, avgLab="US",
                 unitLabX=1995.125, labSet="tb", labK=1,
                 initialZero=TRUE, CI=FALSE,
                 limitsFD=limits, limitsRR=limits, xat=xat,
                 col=col, dots=TRUE, RRasPD=TRUE,
                 plottitles=plottitles,
                 xtitles=xtitles, ytitles=ytitles)

# Plot the results by state with weighted mean
cigUnitLineplots(FD=simFD.3g.scen1, RR=simRR.3g.scen1,
                 periods=periods, units=stateList,
                 showUnits=FALSE, showMean=TRUE,
                 avg=avg, avgLab="US",
                 unitLabX=1995.125, labSet="tb", labK=1,
                 initialZero=TRUE, CI=TRUE,
                 limitsFD=limits, limitsRR=limits, xat=xat,
                 col=col, dots=TRUE, RRasPD=TRUE,
                 plottitles=plottitles,
                 xtitles=xtitles, ytitles=ytitles)

# Plot the representative state results
cigUnitLineplots(FD=simFD.3g.scen1, RR=simRR.3g.scen1,
                 compFD=sim.fd3g, compRR=sim.rr3g,
                 periods=periods, units=stateList,
                 showUnits=FALSE, showMean=FALSE,
                 avg=avg, avgLab="Wgt", compLab="Rep",
                 unitLabX=1995.15, labSet="tb", labK=1,
                 initialZero=TRUE, CI=TRUE, compCI=TRUE,
                 limitsFD=limits, limitsRR=limits, xat=xat,
                 col=col, dots=TRUE, RRasPD=TRUE,
                 plottitles=plottitles,
                 xtitles=xtitles, ytitles=ytitles)

# Plot the comparison of aggregate results and representative results
cigUnitLineplots(FD=simFD.3g.scen1, RR=simRR.3g.scen1,
                 compFD=sim.fd3g, compRR=sim.rr3g,
                 periods=periods, units=stateList,
                 showUnits=FALSE, showMean=TRUE,
                 avg=avg, avgLab="Wgt", compLab="Rep",
                 unitLabX=1995.15, labSet="tb", labK=1,
                 initialZero=TRUE, CI=TRUE, compCI=TRUE,
                 limitsFD=limits, limitsRR=limits, xat=xat,
                 col=col, dots=TRUE, RRasPD=TRUE,
                 plottitles=plottitles,
                 xtitles=xtitles, ytitles=ytitles)

# Plot dotplot of results for final period by state: all states
cigUnitDotplot(x=simFD.3g.scen1$units,
               units=stateFull, period=3,
               showRank=TRUE,
               vertmark=0,
               limits=c(-62, 5),
               labspace=0.35,
               at=seq(-60,0,10),
               col=col[1],
               xtitle="Change, Packs pc 3 years after +$.60 tax",
               toptitle="Change, Packs pc 3 years after +$.60 tax")

# Plot dotplot of results for final period by state: top half
cigUnitDotplot(x=simFD.3g.scen1$units,
               units=stateFull, period=3,
               showRank=TRUE,
               showSelected=1:24,
               vertmark=0,
               limits=c(-62, 5),
               labspace=0.35,
               at=seq(-60,0,10),
               col=col[1],
               xtitle="Change, Packs pc 3 years after +$.60 tax",
               toptitle="Change, Packs pc 3 years after +$.60 tax")

# Plot dotplot of results for final period by state: bottomhalf
cigUnitDotplot(x=simFD.3g.scen1$units,
               units=stateFull, period=3,
               showRank=TRUE,
               showSelected=25:48,
               vertmark=0,
               limits=c(-62, 5),
               labspace=0.35,
               at=seq(-60,0,10),
               col=col[1],
               xtitle="Change, Packs pc 3 years after +$.60 tax",
               toptitle="Change, Packs pc 3 years after +$.60 tax")