# R script used in https://remcat.hatenadiary.jp/entry/20211009/maikinold
# URL: http://tsigeto.info/maikin/maikin-monthly3.r.txt
# 2021/09/26 - 2021/10/05
# Created by TANAKA Sigeto <http://tsigeto.info/maikin/>

datafile.new <- "maikin-monthly.dat"
datafile.old <- "maikin-monthly-juu.dat"
datafile.kyu201801 <- "maikin201801kyu.dat"

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 )

read.data.long <- function( filename ) {
	d <- read.delim( filename, header=T )

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

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

	d1 <- d
	d2 <- d
	d1$worker <- d1$e0
	d2$worker <- d2$e1
	d2$yyyymm <- d2$yyyymm + 0.5
	r <- rbind( d1, d2 )
	r[ order(r$yyyymm) , ]
}
list.bysize <- function( data , reset=c() ){
	r <- list(
		size5 = subset( data,   5 == size ),
		size30= subset( data,  30 == size ),
		size100=subset( data, 100 == size ),
		size500=subset( data, 500 == size ),
		size1000=subset(data,1000 == size )
	)
	lapply(
		r ,
		function(d){
			n <- nrow( d )
			d$worker.prev <- c( NA, d [ -n, "worker" ] )
			d$worker.next <- c(     d [ -1, "worker" ] , NA )
			d$worker.inc <- log( d$worker / d$worker.prev )
			d$worker.diff<- d$worker - d$worker.prev 
			d [ d$yyyymm %in% reset , "worker.inc" ] <- NA
			d [ d$yyyymm %in% reset , "worker.diff"] <- NA
			d
		}
	)
}

# Renewed data
x.long <- read.data.long( datafile.new )
x.bysize <- subset( x.long, industry=="TL" & 0<size )
data.new <- list.bysize( x.bysize, reset=reset.date1 )

# Old data
y.long <- read.data.long( datafile.old )
y.bysize <- subset( y.long, industry=="TL" & 0<size )
data.old <- list.bysize ( y.bysize, reset=reset.date1 )
start.old.mon <- data.old[[1]] [ 1 , "yyyymm" ]
end.old.mon <- data.old[[1]] [ nrow(data.old[[1]]) , "yyyymm" ]
start.old.num <- which( data.new[[1]]$yyyymm==start.old.mon )
skip.old <- start.old.num -1 
remain.old <- nrow( data.new[[1]] ) - nrow( data.old[[1]] ) - skip.old

# Workers at 2012-01
w201201 <- cbind(
	subset( x.bysize, yyyymm==201201 , select=c("yyyymm", "size", "worker") ) ,
	subset( y.bysize, yyyymm==201201 , select="worker" )
)
w201201 <- cbind( w201201, w201201[,3] - w201201[,4] )
colnames(w201201) <- c( "month", "size", "worker.new", "worker.old", "difference" ) 
w201201


# Gaps
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( datafile.kyu201801, 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 <- 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), 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.new[["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.new[["size5"]]$yyyymm

# Trend of the number of workers by establishment size
result <- lapply(
	data.new , 
	function(d){
		size <- d[1,"size"]
		r <- cbind(
			d$yyyymm,
			d$worker,
			worker.pop.yyyymm[ , as.character(size) ]
		)
		colnames(r) <- c( "yyyymm", "worker.new", "census" )
		rownames(r) <- rownames(d)
		data.frame(r)
	}
)
worker.old <- sapply(
	data.old , 
	function(d){
		r <- d$worker
		names(r) <- d$yyyymm
		r
	}
)
for( s in names(result) ) {
	result[[s]]$worker.old <- c( rep(NA, skip.old) , worker.old[, s] , rep(NA, remain.old) )
}

# Select record and field
select.f <- function( data, mod, field ){
	d <- subset( data, yyyymm %% 1 == mod )
	r <- d[ , field]
	names(r) <- floor( d$yyyymm )
	r
}

# Distinguish two factors
worker.cum <- function( v , start ) {
	v[ is.na(v) ] <- 0
	v[ 1 ] <- 0
	start * exp( cumsum(v) )
}

cum5.new <- sapply(
	data.new ,
	function(d) { 
		size <-  d[1, "size"]
		start <- d[1, "worker"]
		worker.cum( select.f( d, 0.5 , "worker.inc" ) , start )
	}
)
cum0.new <- sapply(
	data.new ,
	function(d) { 
		size <-  d[1, "size"]
		start <- d[1, "worker"]
		worker.cum( select.f( d, 0 , "worker.inc" ) , start )
	}
)

# Workers at 2007-10
w200710 <- rbind( cum0.new[ "200710", ], cum5.new[ "200710", ] )
rownames(w200710)<- c( "e1.e0", "e0.e1" )
colnames(w200710)<- c( 5, 30, 100, 500, 1000)

cum5.old <- sapply(
	data.old ,
	function(d) { 
		size <-  d[1, "size"]
		start <- w200710[ "e0.e1",  as.character(size) ]
		worker.cum( select.f( d, 0.5 , "worker.inc" ) , start )
	}
)
cum0.old <- sapply(
	data.old ,
	function(d) { 
		size <-  d[1, "size"]
		start <- w200710[ "e1.e0",  as.character(size) ]
		worker.cum( select.f( d, 0 , "worker.inc" ) , start )
	}
)

# Census
temp <- sapply (
	as.numeric( rownames(cum0.new) ) ,
	function(i) {
		j <- as.character(i)
		k <- as.character(i-0.5)
		w <- rownames(worker.pop)
		if (     j %in% w ) { r <- worker.pop[ j, ] }
		else if( k %in% w ) { r <- worker.pop[ k, ] }
		else { r <- rep( NA, ncol(worker.pop) ) }
		r
	}	
)
worker.pop.yyyymm2 <- t(temp)
colnames(worker.pop.yyyymm2) <- colnames(worker.pop)
rownames(worker.pop.yyyymm2) <- rownames(cum0.new)

worker.predicted <- list()
for( i in 1:ncol(cum0.new) ) {
	s <- colnames(cum0.new)[i]
	worker.predicted[[s]] <- data.frame( yyyymm=rownames(cum0.new) )
	worker.predicted[[s]]$e1.e0.new <- cum0.new[ , i] 
	worker.predicted[[s]]$e0.e1.new <- cum5.new[ , i] 
	worker.predicted[[s]]$census <- worker.pop.yyyymm2[ , i ]
	worker.predicted[[s]]$e1.e0.old <- c( rep(NA,skip.old/2), cum0.old[,i] , rep(NA,remain.old/2) )
	worker.predicted[[s]]$e0.e1.old <- c( rep(NA,skip.old/2), cum5.old[,i] , rep(NA,remain.old/2) )
}

inc.e0.e1.new <- sapply( data.new , select.f , 0.5 , "worker.inc")
inc.e0.e1.old <- sapply( data.old , select.f , 0.5 , "worker.inc" )
inc.e0.e1 <- list()
for( i in 1:ncol(inc.e0.e1.new) ) {
	s <- colnames(inc.e0.e1.new)[i]
	inc.e0.e1[[s]] <- data.frame( yyyymm=rownames(inc.e0.e1.new) )
	inc.e0.e1[[s]]$new <- inc.e0.e1.new[ , i ]
	inc.e0.e1[[s]]$old <- c( rep(NA,skip.old/2), inc.e0.e1.old[,i] , rep(NA,remain.old/2) )
}

inc.e1.e0.new <- sapply( data.new , select.f , 0   , "worker.inc")
inc.e1.e0.old <- sapply( data.old , select.f , 0   , "worker.inc" )
inc.e1.e0 <- list()
for( i in 1:ncol(inc.e1.e0.new) ) {
	s <- colnames(inc.e1.e0.new)[i]
	inc.e1.e0[[s]] <- data.frame( yyyymm=rownames(inc.e1.e0.new) )
	inc.e1.e0[[s]]$new <- inc.e1.e0.new[ , i ]
	inc.e1.e0[[s]]$old <- c( rep(NA,skip.old/2), inc.e1.e0.old[,i] , rep(NA,remain.old/2) )
}