Duration for Letta government

letta.R

## ############################
## LOAD THE NECESSARY LIBRARIES
## ############################
library(foreign)
library(survival)

## ######################
## READ DATA AND DO DATES
## ######################

## Assumes you have the data from http://www.erdda.se/
dat<-read.dta("CII-2007Final.dta")

## start date
dat$v004x<-as.character(dat$v004x)
dat$v004x<-paste("19",dat$v004x,sep="")
dat$startdate<-as.Date(dat$v004x,format="%Y%m%d")

## start date
dat$v005x<-as.character(dat$v005x)
out_of_sample<-which(dat$v005x=="88888")
dat$v005x<-paste("19",dat$v005x,sep="")
dat$v005x[out_of_sample]<-"19990101"
dat$enddate<-as.Date(dat$v005x,format="%Y%m%d")

## #############################
## CALCULATE DURATION AND STATUS
## #############################
dat$durat<-as.numeric(dat$enddate-dat$startdate)
dat$durat[dat$durat<0]<-NA
dat$status<-1
dat$status[out_of_sample]<-0 ## out of sample
dat$status[which(dat$v217y7==1)]<-0 ## technical termination

## OR, alternately
dat$status<-1
dat$status[out_of_sample]<-0 ## out of sample
dat$status[which(dat$v217y==0)]<-0 ## technical termination

## OR, alternately (preferred)
dat$status<-1
dat$status[out_of_sample]<-0 ## out of sample
dat$status[pmax(dat$v217y2,dat$v217y32)==1]<-0 ## technical termination

## ############
## DO COX MODEL
## ############

## Note that:
## breslow for comparability with stata/spss
## FRG is the reference category for countries


my.cox.country_dummies<-coxph(Surv(durat,status)~v001y+v002y+v003y+v004y+v005y+v007y+v008y+v009y+v010y+v011y+v012y+v013y+v014y+v015y+v016y+v017y,data=dat,method="breslow")
my.cox.decade<-update(my.cox.country_dummies,.~.+v031y2 + v033y2 + v034y2 + v035y2 + v036y2)
summary(my.cox.decade)

my.cox.structure<-coxph(Surv(durat,status)~v056y+v050y+v051y+v057y+v052y+v049y+v041y+v043y,data=dat,method="breslow")
my.cox.bestfit<-coxph(Surv(durat,status)~v056y+v057y+v052y+v049y+v041y+v043y+v091y2+v092y+v083y+v126y+v125y+v130y+v134y+v124y2+v137y+v205y,data=dat,method="breslow")

## recode into human-readable names
names(dat)<-car::recode(names(dat),
	"'v056y'='majcab';'v057y'='mwin';'v052y'='maxbpp';'v049y'='coalition';
	'v041y'='maxdurat';'v043y'='enpp';'v091y2'='minconnected';'v092y'='conservative';'v083y'='issdims';'v126y'='pparl';
	'v125y'='oppinf';'v130y'='cabunam';'v134y'='pmdiss';'v124y2'='bicam';'v137y'='semipres';'v205y'='inflation'")

## This taken from the 
my.cox.bestfit<-coxph(Surv(durat,status)~
	majcab+mwin+maxbpp+coalition+maxdurat+enpp+minconnected+conservative+issdims+pparl+oppinf+cabunam+pmdiss+bicam+semipres+v010y,
	data=dat,method="breslow")

my.weibull.bestfit<-survreg(formula(my.cox.bestfit),dist="weibull",data=dat)
my.weibull.bestfit<-update(my.weibull.bestfit,.~.+v017y)

## Now generate predictions
seat.shares <- c(109,7,1,6,99,17,1,54,19,1,1) ## senato
seat.shares<-c(297,37,6,5,98,18,9,109,39,8,2,1,1)  ## camera?

library(Zelig)

z.out<-zelig(formula(my.cox.bestfit),data=dat,model="coxph")


x.out1<-setx(z.out,
	majcab=1, ## the coalition has a majority
	mwin=0, ## the coalition is minimal winning
	maxbpp=1, ## Banzhaf() = {2548, 1548, 1548, 436, 328, 244, 160, 160, 56, 56, 56, 56}
	coalition=1, ## the coalition is a coalition
	maxdurat=(5*365 - 63), ## through inspecting dat$maxdurat[which(dat$v017y==1)]
	enpp= 1/sum((seat.shares/sum(seat.shares))^2), ##
	minconnected=1, ##
	conservative=0, ## i.e., a conservative party has a majority of cabinet seats: see 97fn6 of the book
	issdims=mean(dat$issdims[which(dat$v010y==1)],na.rm=T), ## through inspecting mean(dat$issdims[which(dat$v010y==1)])
	pparl=1, ## investiture vote required
	oppinf=mean(dat$oppinf[which(dat$v010y==1)]), ## dat$oppinf[which(dat$v017y==1)]
	cabunam=1, ## unanimity required in cabinet
	pmdiss=0, ## PM has dismissal powers
	bicam=1, ## There's the Senato
	semipres=0,  ## no semipres
	v010y=1) ## 



s.out1<-sim(z.out,x.out1)

## 

foo <- summary(s.out1)
x <- as.numeric(rownames(foo$qi.stats$survival))
y <- foo$qi.stats$survival[,1]
ymax <-foo$qi.stats$survival[,4]
ymin <-foo$qi.stats$survival[,3]  
plot.df <- data.frame(x=x,y=y,ymax=ymax,ymin=ymin)
the.day <- x[which.min(abs(y-.5))]
plot.df$the.day <- the.day

library(ggplot2)
png(file="letta_durat.png",width=700,height=480)
ggplot(plot.df,aes(x=x,y=y,ymax=ymax,ymin=ymin)) + 
	geom_smooth(stat="identity") + theme_bw() +
	scale_x_continuous("Days") + 
	scale_y_continuous("Prob. of surviving") +
	geom_vline(xintercept = the.day,lty=2) + geom_hline(yintercept = 0.5,lty =2)
dev.off()