# R script used in https://remcat.hatenadiary.jp/entry/20210911/gap
# URL: http://tsigeto.info/maikin/maikin-monthly.r.txt
# 2021/08/16 - 2021/09/11
# Created by TANAKA Sigeto <http://tsigeto.info/maikin/>

x <- read.delim( "maikin-monthly.dat", header=T )

# Sort by date and establishment size
x <- x[ order(x$size) , ]
x <- x[ order(x$yyyymm) , ]

x$is <- factor( paste( x$industry, x$size, sep="." ) )
x$year <- round( x$yyyymm / 100 )
x$month<- floor( x$yyyymm %% 100 )
x$worker2 <- ( x$e0 + x$e1 ) /2

x1 <- x
x2 <- x
x1$worker <- x1$e0
x2$worker <- x2$e1
x2$yyyymm <- x2$yyyymm + 0.5
temp <- rbind( x1, x2 )
x.long <- temp[ order(temp$yyyymm) , ]

x.bysize <- subset( x.long, industry=="TL" & 0<size )
data <- list(
	size5 = x.bysize[  5 == x.bysize$size , ],
	size30= x.bysize[ 30 == x.bysize$size , ],
	size100=x.bysize[100 == x.bysize$size , ],
	size500=x.bysize[500 == x.bysize$size , ],
	size1000=x.bysize[1000==x.bysize$size , ]
)

subset( data[["size5"]], 201901<=yyyymm, select=c( "yyyymm", "worker" ) )

census.date <- c( 200610, 200907, 201407, 201605.5, 201905.5 )

reset.date0 <- c(  200812.5, 201112.5, 201712.5 )
reset.date1 <- c(  200901  , 201201  , 201801   )
checkpoint <- c( census.date, reset.date0, reset.date1, 202105, 202105.5 )
x.cp <- subset( x.bysize, yyyymm %in% checkpoint )
size.cp <- sapply( split( x.cp , x.cp$yyyymm ), function(d){ d[ , "size" ] } ) 
worker.cp<-sapply( split( x.cp , x.cp$yyyymm ), function(d){ d[ , "worker" ] } )
rownames(worker.cp) <- size.cp[,1]

gap2009 <- worker.cp[,"200901"]/ worker.cp[,"200812.5"]
worker.pop2006 <- worker.cp[,"200610"] * gap2009

gap2012 <- worker.cp[,"201201"]/ worker.cp[,"201112.5"]
worker.pop2009 <- worker.cp[,"200907"] * gap2012

kyu2018 <- read.delim( "maikin201801kyu.dat", header=T )
kyu2018 <- kyu2018[ order(kyu2018$size),  ]
kyu2018.tl0 <- subset( kyu2018, industry=="TL" & 0<size )
gap2018 <- worker.cp[,"201801"] / kyu2018.tl0$e0
worker.pop2014 <- worker.cp[,"201407"] * gap2018

# From https://www.mhlw.go.jp/content/10700000/000823050.pdf
gap2019 <- c( 0.888, 1.092, 0.985, 0.813, 0.947 )
names( gap2019 ) <- c( "5", "30",  "100",  "500",  "1000"  )
worker.pop2019 <- worker.cp[,"201905.5"] * gap2019

census2016.rev <- rbind (
	c( 0.878, 1.161, 0.960, 0.843, 0.971 ) ,
	c( 0.880, 1.171, 0.960, 0.844, 0.973 ) ,
	c( 0.877, 1.162, 0.959, 0.843, 0.977 ) ,
	c( 0.877, 1.161, 0.960, 0.842, 0.970 ) ,
	c( 0.877, 1.161, 0.960, 0.842, 0.970 )
)
gap2016 <- exp( apply( log(census2016.rev), 2, mean ) )
names( gap2016 ) <- c( "5", "30",  "100",  "500",  "1000"  )
worker.pop2016 <- worker.cp[,"201605.5"] * gap2016

worker.pop <- rbind(
	worker.pop2006,
	worker.pop2009,
	worker.pop2014,
	worker.pop2016,
	worker.pop2019
)
rownames(worker.pop) <- census.date

temp <- sapply (
	data[["size5"]]$yyyymm ,
	function(i) {
		if( i %in% rownames(worker.pop) ) { r <- worker.pop[ as.character(i), ] }
		else { r <- rep( NA, ncol(worker.pop) ) }
		r
	}	
)
worker.pop.yyyymm <- t(temp)
colnames(worker.pop.yyyymm) <- colnames(worker.pop)
rownames(worker.pop.yyyymm) <- data[["size5"]]$yyyymm

result <- lapply(
	data , 
	function(d){
		size <- d[1,"size"]
		r <- cbind(
			d$yyyymm,
			d$worker,
			worker.pop.yyyymm[ , as.character(size) ]
		)
		colnames(r) <- c( "yyyymm", "worker", "census" )
		rownames(r) <- rownames(d)
		data.frame(r)
	}
)

noinc <- lapply(
	data, 
	function(d) {
		size <- d[1,"size"]
		d$census <- worker.pop.yyyymm[ , as.character(size) ]
		d$noinc <- NA
		d[1, "noinc"] <- d[1, "worker"]
		d[ d$yyyymm==200901, "noinc" ] <- worker.pop[ "200610"  , as.character(size) ]
		d[ d$yyyymm==201201, "noinc" ] <- worker.pop[ "200907"  , as.character(size) ]
		d[ d$yyyymm==201801, "noinc" ] <- worker.pop[ "201407"  , as.character(size) ]
		for( i in 2:nrow(d) ) {
			if( is.na( d[ i, "noinc"] ) ) { d[ i, "noinc"] <- d[ i-1, "noinc"] }
		}
		d
	}
)

result2 <- lapply(
	noinc , 
	function(d){
		d$worker <- d$noinc
		d[ , c( "yyyymm", "worker", "census" ) ]
	}
)

worker.census.gap <- function(d){
	c <- subset( d, !is.na(census) )
	c$gap <- c$worker / c$census
	c
}
result.gap <- lapply( result, worker.census.gap )
result2.gap<-lapply( result2, worker.census.gap )

gap.table <- data.frame()
for( s in names(result.gap) ) {
	r <- result.gap[[s]]
	r$is <- s
	r$noinc <- result2.gap[[s]]$worker
	r$gap2 <- result2.gap[[s]]$gap
	r$gap.rate <- r$gap / r$gap2
	r$diff1 <- r$worker - r$census
	r$diff2 <- r$noinc - r$census
	r$diff.diff <- abs( r$diff1 ) - abs( r$diff2 )
	gap.table <- rbind( gap.table, r )
}

temp <- split( gap.table, gap.table$is )
diff.diff <- sapply( temp, function(d){ d[, "diff.diff"] } )
rownames(diff.diff) <- result.gap[["size5"]]$yyyymm
diff.diff <- diff.diff[ , names(result.gap) ]
round( addmargins(diff.diff) )