head(aa)
# Last year number of points, goals scored and conceived in the league
SeasonID <- unique(dataset$Season)
for(s in 2:length(SeasonID)){
df <- dataset %>% filter(Season==SeasonID[s-1])
leaguetable <- maketable(df)
merge(data.frame("team"=teams),leaguetable,by="team",all.x=T)
# home
aa <- merge(dataset[dataset$Season==SeasonID[s],],leaguetable[,c("team","GP","gf","ga","Pts")],
by.x="home",by.y="team",all.x=T)
aa <- aa[order(aa[,"ID"]),]
dataset[dataset$Season==SeasonID[s],"hlastyrpts"] <- aa[,"Pts"]
dataset[dataset$Season==SeasonID[s],"hlasyrgfpergame"] <- aa[,"gf"]/aa[,"GP"]
dataset[dataset$Season==SeasonID[s],"hlasyrgapergame"] <- aa[,"ga"]/aa[,"GP"]
# visitor
aa <- merge(dataset[dataset$Season==SeasonID[s],],leaguetable[,c("team","GP","gf","ga","Pts")],
by.x="visitor",by.y="team",all.x=T)
aa <- aa[order(aa[,"ID"]),]
dataset[dataset$Season==SeasonID[s],"vlastyrpts"] <- aa[,"Pts"]
dataset[dataset$Season==SeasonID[s],"vlasyrgfpergame"] <- aa[,"gf"]/aa[,"GP"]
dataset[dataset$Season==SeasonID[s],"vlasyrgapergame"] <- aa[,"ga"]/aa[,"GP"]
}
dataset[is.na(dataset[,"hlastyrpts"]),"hlastyrpts"] <- 0
dataset[is.na(dataset[,"vlastyrpts"]),"vlastyrpts"] <- 0
head(dataset)
i=1000
datatset[i,]
dataset[i,]
newdata <- filter(dataset,home==home[i,"ID"] && visitor==visitor[i,"ID"])
home[i,"ID"]
dataset[i,"home"]
dataset[i,"visitor"]
newdata <- filter(dataset,home==dataset[i,"home"] & visitor==dataset[i,"visitor"])
newdata
newdata <- filter(dataset,home==dataset[i,"home"] & visitor==dataset[i,"visitor"] & ID < dataset[i,"ID"])
newdata
newdata <- filter(dataset,home==dataset[i,"home"] & visitor==dataset[i,"visitor"] & ID < 600)
newdata
# Previous game between the two teams
GetLstGameScore <- function(i,dataset){
newdata <- filter(dataset,home==dataset[i,"home"] & visitor==dataset[i,"visitor"] & ID < dataset[i,"ID"])
hScore <- newdata[nrow(newdata),"hgoal"]
vScore <- newdata[nrow(newdata),"vgoal"]
return(list(hScore=hScore,vScore=vScore))
}
GetLstGameScore(1000,dataset)
GetLstGameScore(10,dataset)
is.na(GetLstGameScore(10,dataset))
is.na(GetLstGameScore(10,dataset)$hScore)
is.integer(GetLstGameScore(10,dataset)$hScore)
GetLstGameScore(10,dataset)$hScore < 0
GetLstGameScore(10,dataset)$hScore==integer(0)
GetLstGameScore <- function(i,dataset){
newdata <- filter(dataset,home==dataset[i,"home"] & visitor==dataset[i,"visitor"] & ID < dataset[i,"ID"])
if(nrow(newdata)==0){
hScore <- NA
vScore <- NA
} else {
hScore <- newdata[nrow(newdata),"hgoal"]
vScore <- newdata[nrow(newdata),"vgoal"]
}
return(list(hScore=hScore,vScore=vScore))
}
GetLstGameScore(1000,dataset)
GetLstGameScore(10,dataset)
# Previous game between the two teams
GetLstGameScore <- function(i,dataset){
newdata <- filter(dataset,home==dataset[i,"home"] & visitor==dataset[i,"visitor"] & ID < dataset[i,"ID"])
if(nrow(newdata)==0){
Score <- c(NA,NA)
} else {
Score <- c(newdata[nrow(newdata),"hgoal"],newdata[nrow(newdata),"vgoal"])
}
return(Score)
}
GetLstGameScore(1000,dataset)
GetLstGameScore(10,dataset)
aa <- mapply(GetLstGameScore,1:nrow(dataset),dataset)
? apply
? mapply
aa <- mapply(GetLstGameScore,1:2,dataset)
mapply(GetLstGameScore,1:2,dataset)
mapply(GetLstGameScore,1,dataset)
GetLstGameScore(1,dataset)
for(i in 1:nrow(dataset)){
print(i)
dataset[,c("LstGamehgoal","LstGamevScore")] <- GetLstGameScore(i,dataset)
}
head(dataset)
dataset[995:1000,]
GetLstGameScore(995,dataset)
GetLstGameScore(996,dataset)
GetLstGameScore(997,dataset)
for(i in 1:nrow(dataset)){
print(i)
Score <- GetLstGameScore(i,dataset)
dataset[,"LstGamehgoal"] <- Score[1]
dataset[,"LstGamevgoal"] <- Score[2]
}
head(dataset)
1000%%3
10%%3
9%%3
for(i in 1:nrow(dataset)){
if(i%%500==0) print(i)
Score <- GetLstGameScore(i,dataset)
dataset[,"LstGamehgoal"] <- Score[1]
dataset[,"LstGamevgoal"] <- Score[2]
}
tail(dataset)
GetLstGameScore(24644,dataset)
GetLstGameScore(24643,dataset)
GetLstGameScore(24642,dataset)
dataset[1000,c("LstGamehgoal","LstGamevgoal")]
dataset[1000:1020,c("LstGamehgoal","LstGamevgoal")]
Score
i
GetLstGameScore(i,dataset)
rm(list=ls())
### Set working directory
setwd("R:/Simulations/BEAST")
rm(list=ls())
### Set working directory
setwd("Volumes/USB_KEY:/Simulations/BEAST")
### Load packages
library("Synth")
install.packges("Synth")
install.packages("Synth")
library("Synth")
setwd("R:/Simulations/BEAST")
setwd("/Users/jeremylhour/Documents/R/BEAST")
data <- data.frame(t(read.table("E:/BEAST/datasets/MLAB_data.txt")))
#data <- data.frame(t(read.table("Z:/Simulations/MLAB_data.txt")))
Names <- c("State_ID","Income","Retail price", "Young", "BeerCons","Smoking1988", "Smoking1980","Smoking1975",
mapply(function(x) paste("SmokingCons",x,sep=""),1970:2000))
colnames(data) <- Names
data[,"Treated"] <- as.numeric(data[,"State_ID"]==3) #California is state with ID=3
CaliSmoke <- ts(unlist(data[data[,"Treated"]==1, mapply(function(x) paste("SmokingCons",x,sep=""),1970:2000)]),
start=c(1970), freq=1)
data <- data.frame(t(read.table("/datasets/MLAB_data.txt")))
#data <- data.frame(t(read.table("Z:/Simulations/MLAB_data.txt")))
data <- data.frame(t(read.table("/datasets/MLAB_data.txt")))
#data <- data.frame(t(read.table("Z:/Simulations/MLAB_data.txt")))
setwd("/Users/jeremylhour/Documents/R/BEAST")
data <- data.frame(t(read.table("/datasets/MLAB_data.txt")))
#data <- data.frame(t(read.table("Z:/Simulations/MLAB_data.txt")))
data <- data.frame(t(read.table("/datasets/MLAB_data.txt")))
data <- data.frame(t(read.table("datasets/MLAB_data.txt")))
data
Names <- c("State_ID","Income","Retail price", "Young", "BeerCons","Smoking1988", "Smoking1980","Smoking1975",
mapply(function(x) paste("SmokingCons",x,sep=""),1970:2000))
colnames(data) <- Names
data[,"Treated"] <- as.numeric(data[,"State_ID"]==3) #California is state with ID=3
CaliSmoke <- ts(unlist(data[data[,"Treated"]==1, mapply(function(x) paste("SmokingCons",x,sep=""),1970:2000)]),
start=c(1970), freq=1)
### Figure 1: California cigarettes consumption
plot(CaliSmoke, col="steelblue", lwd=2,
xlab="", ylab="Packs per capita",
ylim=c(0,150),
main="California cigarette pack consumption")
abline(v=1988, col="firebrick")
# Pre-treatment outcome variables to optimize the fit
Z = NULL
for(t in c(1970:1975, 1980, 1988)){
varname <- paste("SmokingCons",t,sep="")
Z <- cbind(Z,data[,varname])
}
# Find the weights...
ADH_SC <- synth(data.prep.obj = NULL,
X1 = as.matrix(X[d==1,-1]), X0 = t(as.matrix(X[d==0,-1])),
Z1 = as.matrix(Z[d==1,]), Z0 = t(as.matrix(Z[d==0,])),
custom.v = NULL,
optimxmethod = "Nelder-Mead",
genoud = FALSE, quadopt = "ipop",
Margin.ipop = 5e-04,
Sigf.ipop = 5,
Bound.ipop = 10,
verbose = FALSE)
### Figure 1: California cigarettes consumption
plot(CaliSmoke, col="steelblue", lwd=2,
xlab="", ylab="Packs per capita",
ylim=c(0,150),
main="California cigarette pack consumption")
abline(v=1988, col="firebrick")
### Treatment effect
X <- cbind(rep(1,nrow(data)),
data[,c("Income","Retail price", "Young", "BeerCons"
, "SmokingCons1970", "SmokingCons1971", "SmokingCons1972", "SmokingCons1973", "SmokingCons1974", "SmokingCons1975"
#, "SmokingCons1976", "SmokingCons1977", "SmokingCons1978", "SmokingCons1979"
, "SmokingCons1980"
#, "SmokingCons1981", "SmokingCons1982", "SmokingCons1983", "SmokingCons1984"
#, "SmokingCons1985", "SmokingCons1986", "SmokingCons1987"
, "SmokingCons1988"
)])
d <- data[,"Treated"]
# Pre-treatment outcome variables to optimize the fit
Z = NULL
for(t in c(1970:1975, 1980, 1988)){
varname <- paste("SmokingCons",t,sep="")
Z <- cbind(Z,data[,varname])
}
# Find the weights...
ADH_SC <- synth(data.prep.obj = NULL,
X1 = as.matrix(X[d==1,-1]), X0 = t(as.matrix(X[d==0,-1])),
Z1 = as.matrix(Z[d==1,]), Z0 = t(as.matrix(Z[d==0,])),
custom.v = NULL,
optimxmethod = "Nelder-Mead",
genoud = FALSE, quadopt = "ipop",
Margin.ipop = 5e-04,
Sigf.ipop = 5,
Bound.ipop = 10,
verbose = FALSE)
rm(list=ls())
### Load packages
library("Synth")
data <- data.frame(t(read.table("datasets/MLAB_data.txt")))
Names <- c("State_ID","Income","Retail price", "Young", "BeerCons","Smoking1988", "Smoking1980","Smoking1975",
mapply(function(x) paste("SmokingCons",x,sep=""),1970:2000))
colnames(data) <- Names
data[,"Treated"] <- as.numeric(data[,"State_ID"]==3) #California is state with ID=3
CaliSmoke <- ts(unlist(data[data[,"Treated"]==1, mapply(function(x) paste("SmokingCons",x,sep=""),1970:2000)]),
start=c(1970), freq=1)
### Figure 1: California cigarettes consumption
plot(CaliSmoke, col="steelblue", lwd=2,
xlab="", ylab="Packs per capita",
ylim=c(0,150),
main="California cigarette pack consumption")
abline(v=1988, col="firebrick")
### Treatment effect
X <- cbind(rep(1,nrow(data)),
data[,c("Income","Retail price", "Young", "BeerCons"
, "SmokingCons1970", "SmokingCons1971", "SmokingCons1972", "SmokingCons1973", "SmokingCons1974", "SmokingCons1975"
#, "SmokingCons1976", "SmokingCons1977", "SmokingCons1978", "SmokingCons1979"
, "SmokingCons1980"
#, "SmokingCons1981", "SmokingCons1982", "SmokingCons1983", "SmokingCons1984"
#, "SmokingCons1985", "SmokingCons1986", "SmokingCons1987"
, "SmokingCons1988"
)])
d <- data[,"Treated"]
### Setting
d <- as.matrix(d)
X <- as.matrix(X)
n <- nrow(X)
p <- ncol(X)
# Pre-treatment outcome variables to optimize the fit
Z = NULL
for(t in c(1970:1975, 1980, 1988)){
varname <- paste("SmokingCons",t,sep="")
Z <- cbind(Z,data[,varname])
}
# Find the weights...
ADH_SC <- synth(data.prep.obj = NULL,
X1 = as.matrix(X[d==1,-1]), X0 = t(as.matrix(X[d==0,-1])),
Z1 = as.matrix(Z[d==1,]), Z0 = t(as.matrix(Z[d==0,])),
custom.v = NULL,
optimxmethod = "Nelder-Mead",
genoud = FALSE, quadopt = "ipop",
Margin.ipop = 5e-04,
Sigf.ipop = 5,
Bound.ipop = 10,
verbose = FALSE)
ADH_SC$solution.w
data
head(data)
i=1
ControlTS <- ts(unlist(data[i, mapply(function(x) paste("SmokingCons",x,sep=""),1970:2000)]),
start=c(1970), freq=1)
ControlTS
plot(ControlTS)
i=2
ControlTS <- ts(unlist(data[i, mapply(function(x) paste("SmokingCons",x,sep=""),1970:2000)]),
start=c(1970), freq=1)
plot(ControlTS)
plot(CaliSmoke)
i=5
ControlTS <- ts(unlist(data[i, mapply(function(x) paste("SmokingCons",x,sep=""),1970:2000)]),
start=c(1970), freq=1)
plot(CaliSmoke)
plot(ControlTS)
lm(CaliSmoke ~ ControlTS)
summary( lm(CaliSmoke ~ ControlTS))
coef(lm(CaliSmoke ~ ControlTS))[2]
a_reg <- vector[, length=1:39]
a_reg <- vector(, length=1:39)
a_reg <- vector(, lenght=1:39)
a_reg <- vector(, lenght=39)
a_reg <- vector(, length=39)
a_reg <- vector(, length=39)
for(i in 1:39){
ControlTS <- ts(unlist(data[i, mapply(function(x) paste("SmokingCons",x,sep=""),1970:2000)]),
start=c(1970), freq=1)
a_reg[i] <- coef(lm(CaliSmoke ~ ControlTS))[2]
}
a_reg
plot(a_reg[-39],ADH_SC$solution.w)
scale(ControlTS)
ControlTS <- scale(ControlTS)
coef(lm(CaliSmoke ~ ControlTS))[2]
coef(lm(scale(CaliSmoke) ~ ControlTS))[2]
corr(scale(CaliSmoke),ControlTS)
cor(scale(CaliSmoke),ControlTS)
a_reg <- vector(, length=39)
for(i in 1:39){
ControlTS <- ts(unlist(data[i, mapply(function(x) paste("SmokingCons",x,sep=""),1970:2000)]),
start=c(1970), freq=1)
a_reg[i] <- cor(CaliSmoke,ControlTS)
}
plot(a_reg[-39],ADH_SC$solution.w)
max(a_reg[-39])
which(a_reg==max(a_reg[-39]))
i=8
ADH_SC$solution.w[8]
ControlTS <- ts(unlist(data[i, mapply(function(x) paste("SmokingCons",x,sep=""),1970:2000)]),
start=c(1970), freq=1)
plot(ControlTS)
plot(cbind(ControlTS,CaliSmoke))
plot(cbind(ControlTS,CaliSmoke), plot.type="single")
### Comparaison des poids avec la correlation des series pre-traitement
a_reg <- vector(, length=39)
for(i in 1:39){
ControlTS <- ts(unlist(data[i, mapply(function(x) paste("SmokingCons",x,sep=""),1970:1988)]),
start=c(1970), freq=1)
a_reg[i] <- cor(CaliSmoke,ControlTS)
}
a_reg <- vector(, length=39)
CaliSmoke <- ts(unlist(data[data[,"Treated"]==1, mapply(function(x) paste("SmokingCons",x,sep=""),1970:1988)]),
start=c(1970), freq=1)
for(i in 1:39){
ControlTS <- ts(unlist(data[i, mapply(function(x) paste("SmokingCons",x,sep=""),1970:1988)]),
start=c(1970), freq=1)
a_reg[i] <- cor(CaliSmoke,ControlTS)
}
plot(a_reg[-39],ADH_SC$solution.w)
max(a_reg)
max(a_reg[-39])
which(a_reg=max(a_reg[-39]))
which(a_reg==max(a_reg[-39]))
i=which(a_reg==max(a_reg[-39]))
ControlTS <- ts(unlist(data[i, mapply(function(x) paste("SmokingCons",x,sep=""),1970:1988)]),
start=c(1970), freq=1)
plot(ControlTS)
plot(cbind(CaliSmoke,ControlTS), plot.type="single")
plot(cbind(scale(CaliSmoke),scale(ControlTS)), plot.type="single")
plot(cbind(scale(CaliSmoke),scale(ControlTS)), plot.type="single", type="line")
scale(CaliSmoke)
scale(ControlTS)
cbind(scale(CaliSmoke),scale(ControlTS))
plot(ts(cbind(scale(CaliSmoke),scale(ControlTS))start=c(1970), freq=1), plot.type="single", type="line")
plot(ts(cbind(scale(CaliSmoke),scale(ControlTS)),start=c(1970), freq=1), plot.type="single", type="line")
CaliSmoke <- ts(unlist(data[data[,"Treated"]==1, mapply(function(x) paste("SmokingCons",x,sep=""),1970:2000)]),
start=c(1970), freq=1)
ControlTS <- ts(unlist(data[i, mapply(function(x) paste("SmokingCons",x,sep=""),1970:2000)]),
start=c(1970), freq=1)
plot(ts(cbind(scale(CaliSmoke),scale(ControlTS)),start=c(1970), freq=1), plot.type="single", type="line")
ControlTS <- ts(unlist(data[i, mapply(function(x) paste("SmokingCons",x,sep=""),1970:1988)]),
start=c(1970), freq=1)
plot(ts(cbind(scale(CaliSmoke),scale(ControlTS)),start=c(1970), freq=1), plot.type="single", type="line")
scale(ControlTS)
scale(ControlTS)$center
center(scale(ControlTS))
m <- mean(scale(ControlTS))
m
m <- mean(ControlTS)
s <- sd(ControlTS)
s
m
(ControlTS-m)/s
plot(ts(cbind(CaliSmoke,mc+sc*(ControlTS-m)/s),start=c(1970), freq=1), plot.type="single", type="line")
sc <- sd(CaliSmoke)
mc <- mean(CaliSmoke)
sc <- sd(CaliSmoke)
mc <- mean(CaliSmoke)
plot(ts(cbind(CaliSmoke,mc+sc*(ControlTS-m)/s),start=c(1970), freq=1), plot.type="single", type="line")
### Comparaison des poids avec la correlation des series pre-traitement
a_reg <- vector(, length=39)
CaliSmoke <- ts(unlist(data[data[,"Treated"]==1, mapply(function(x) paste("SmokingCons",x,sep=""),1970:1988)]),
start=c(1970), freq=1)
for(i in 1:39){
ControlTS <- ts(unlist(data[i, mapply(function(x) paste("SmokingCons",x,sep=""),1970:1988)]),
start=c(1970), freq=1)
a_reg[i] <- cor(CaliSmoke,ControlTS)
}
plot(a_reg[-39],ADH_SC$solution.w)
i=which(a_reg==max(a_reg[-39]))
### Watch pre-treatment fit
ControlTS <- ts(unlist(data[i, mapply(function(x) paste("SmokingCons",x,sep=""),1970:1988)]),
start=c(1970), freq=1)
plot(ts(cbind(scale(CaliSmoke),scale(ControlTS)),start=c(1970), freq=1), plot.type="single", type="line")
m <- mean(ControlTS)
s <- sd(ControlTS)
sc <- sd(CaliSmoke)
mc <- mean(CaliSmoke)
### Watch effect
CaliSmoke <- ts(unlist(data[data[,"Treated"]==1, mapply(function(x) paste("SmokingCons",x,sep=""),1970:2000)]),
start=c(1970), freq=1)
ControlTS <- ts(unlist(data[i, mapply(function(x) paste("SmokingCons",x,sep=""),1970:2000)]),
start=c(1970), freq=1)
plot(ts(cbind(CaliSmoke,mc+sc*(ControlTS-m)/s),start=c(1970), freq=1), plot.type="single", type="line")
i
ADH_SC$solution.w[i]
ADH_SC$solution.w
wei <- ifelse(ADH_SC$solution.w>.1, ADH_SC$solution.w, 0)
wei
i=33
### Watch pre-treatment fit
ControlTS <- ts(unlist(data[i, mapply(function(x) paste("SmokingCons",x,sep=""),1970:1988)]),
start=c(1970), freq=1)
plot(ts(cbind(scale(CaliSmoke),scale(ControlTS)),start=c(1970), freq=1), plot.type="single", type="line")
m <- mean(ControlTS)
s <- sd(ControlTS)
sc <- sd(CaliSmoke)
mc <- mean(CaliSmoke)
CaliSmoke <- ts(unlist(data[data[,"Treated"]==1, mapply(function(x) paste("SmokingCons",x,sep=""),1970:1988)]),
start=c(1970), freq=1)
### Watch pre-treatment fit
ControlTS <- ts(unlist(data[i, mapply(function(x) paste("SmokingCons",x,sep=""),1970:1988)]),
start=c(1970), freq=1)
plot(ts(cbind(scale(CaliSmoke),scale(ControlTS)),start=c(1970), freq=1), plot.type="single", type="line")
m <- mean(ControlTS)
s <- sd(ControlTS)
sc <- sd(CaliSmoke)
mc <- mean(CaliSmoke)
### Watch effect
CaliSmoke <- ts(unlist(data[data[,"Treated"]==1, mapply(function(x) paste("SmokingCons",x,sep=""),1970:2000)]),
start=c(1970), freq=1)
ControlTS <- ts(unlist(data[i, mapply(function(x) paste("SmokingCons",x,sep=""),1970:2000)]),
start=c(1970), freq=1)
plot(ts(cbind(CaliSmoke,mc+sc*(ControlTS-m)/s),start=c(1970), freq=1), plot.type="single", type="line")
NormalSynth <- synth(data.prep.obj = NULL,
X1 = as.matrix(X[d==1,-1]), X0 = t(as.matrix(X[d==0,-1])),
Z1 = as.matrix(X[d==1,-1]), Z0 = t(as.matrix(X[d==0,-1])),
custom.v = rep(1,p-1),
optimxmethod = "Nelder-Mead",
genoud = FALSE, quadopt = "ipop",
Margin.ipop = 5e-04,
Sigf.ipop = 5,
Bound.ipop = 10,
verbose = FALSE)
# Construct counterfactual...
ADH <- vector(length=length(1970:2000))
SC <- vector(length=length(1970:2000))
i <- 0
for(t in 1970:2000){
i=i+1
varname <- paste("SmokingCons",t,sep="")
y <- data[,varname]
ADH[i] <- t(y[d==0])%*%ADH_SC$solution.w
SC[i] <- t(y[d==0])%*%NormalSynth$solution.w
}
### Figure 4: Synthetic Control plot
plotdata <- ts(cbind(CaliSmoke,Immunized,ADH,SC),start=c(1970), freq=1)
### Figure 4: Synthetic Control plot
plotdata <- ts(cbind(CaliSmoke,ADH,SC),start=c(1970), freq=1)
plot(plotdata, plot.type="single",
col=c("steelblue","firebrick","forestgreen","darkgoldenrod1"), lwd=2,
lty=c(1,6,6),xlab="", ylab="Cigarette consumption (Packs per capita)",
ylim=c(35,150))
lim <- par("usr")
rect(1988, lim[3], lim[2], lim[4], col = rgb(0.5,0.5,0.5,1/4))
axis(1) ## add axes back
axis(2)
box()
legend(1971,80,
legend=c("Real California","Synthetic California (weighted)", "Synthetic California (unweighted)"),
col=c("steelblue","firebrick","forestgreen","darkgoldenrod1"), lwd=2,
lty=c(1,6,6))
? synth
plot(plotdata, plot.type="single",
col=c("steelblue","firebrick","forestgreen","darkgoldenrod1"), lwd=2,
lty=c(1,6,6,6),xlab="", ylab="Cigarette consumption (Packs per capita)",
ylim=c(35,150))
lim <- par("usr")
rect(1988, lim[3], lim[2], lim[4], col = rgb(0.5,0.5,0.5,1/4))
axis(1) ## add axes back
axis(2)
box()
legend(1971,80,
legend=c("Real California", "Immunized, Lasso","Synthetic California (weighted)", "Synthetic California (unweighted)"),
col=c("steelblue","firebrick","forestgreen","darkgoldenrod1"), lwd=2,
lty=c(1,6,6,6))
# load data
data(synth.data)
# create matrices from panel data that provide inputs for synth()
dataprep.out<-
dataprep(
foo = synth.data,
predictors = c("X1", "X2", "X3"),
predictors.op = "mean",
dependent = "Y",
unit.variable = "unit.num",
time.variable = "year",
special.predictors = list(
list("Y", 1991, "mean"),
list("Y", 1985, "mean"),
list("Y", 1980, "mean")
),
treatment.identifier = 7,
controls.identifier = c(29, 2, 13, 17, 32, 38),
time.predictors.prior = c(1984:1989),
time.optimize.ssr = c(1984:1990),
unit.names.variable = "name",
time.plot = 1984:1996
)
## run the synth command to identify the weights
## that create the best possible synthetic
## control unit for the treated.
synth.out <- synth(dataprep.out)
## there are two ways to summarize the results
## we can either access the output from synth.out directly
round(synth.out$solution.w,2)
# contains the unit weights or
synth.out$solution.v