# R Code for Construction of Employment Transition Probability Matrices Using Data
# from the 2011 SLID PUMF 
# Created: 12/29/17; Last Modified 04/20/19

# Introduction & Rationale ------------------------------------------------

# NOTE 042019: Modified detection of script location. For this to work, the user
# must set the default working directory under Global Options to the base
# directory of the folder in which the decision model is stored.
# 
# NOTE 111818: The code now detects the path where the script is located, eliminating 
# the need for the user to change this manually upon downloading the script for the 
# first time.
# 
# NOTE 090618: Updated file to draw directly from original PUMF file.

# Initialization ----------------------------------------------------------

# The next line of code clears the variable list, if required.
rm(list=ls())

# Call requisite libraries
library(dplyr); library(readxl); library(data.table); library(here)

Model_Path <- here()
# String for date last modified
LM <- "_042019"
# Filename for output
Filename <- "SLID_2011_Calculations"
# Paths for output
RDA <- paste0(Model_Path, "/R_Objects/")
Data_Path <- paste0(Model_Path, "/Raw_Data/")

# Data Processing ---------------------------------------------------------

# Load the SLID 2011 PUMF
SLID_2011_122917 <- fread(paste0(Data_Path, "slid-75M0010XCB-E-2011pr_F1.csv"))
# This adds a new column that assigns an integer value to age values within a certain 
# range.
SLID_2011_122917 <- mutate(SLID_2011_122917, age_cat = 
                             ifelse(ecage26<25,0,ifelse(ecage26<55,1,2)))
# Similarly, this inserts a column that assigns an integer value to particular levels 
# of educational attainment (based on the value in hleveg18).
SLID_2011_122917 <- mutate(SLID_2011_122917, educ_at = 
                             ifelse(hleveg18>12,99,ifelse(hleveg18>8,1,0)))
# Classify as integers those variables that will be used in the analysis
SLID_2011_122917[c("ecsex99", "age_cat", "educ_at", "ml01v28", "ml11v28", "ml12v28")] <-
  sapply(SLID_2011_122917[c("ecsex99", "age_cat", "educ_at", "ml01v28", 
                            "ml11v28", "ml12v28")], as.integer)

# Use DPLYR to count the number of individuals who were employed, unemployed, or out of
# the labour force in December 2011, grouped by gender, age, education, and labour 
# force status in January 2011.
SLID_2011_1S1 <- SLID_2011_122917 %>% 
  group_by(ecsex99,age_cat,educ_at,ml01v28) %>% 
  summarise("Observations"=n(),"Employed (DEC)"=sum(ml12v28==10),
            "Unemployed (DEC)"=sum(ml12v28==20),"Out of LF (DEC)"=sum(ml12v28==30),
            "N/A (DEC)"=sum(ml12v28==99))
# Now do the same thing, but group by gender, age, education, and labour force status 
# in November 2011.
SLID_2011_2S1 <- SLID_2011_122917 %>% 
  group_by(ecsex99,age_cat,educ_at,ml11v28) %>% 
  summarise("Observations"=n(),"Employed (DEC)"=sum(ml12v28==10),
            "Unemployed (DEC)"=sum(ml12v28==20),"Out of LF (DEC)"=sum(ml12v28==30),
            "N/A (DEC)"=sum(ml12v28==99))

# Here we filter out observations on individuals who were younger than 25 at the 
# beginning of the survey period or who did not know their level of educational
# attainment, as well as those for whom labour force status was described as "not
# applicable" (N=6861) or were out of the labour force in January 2011 or November 
# 2011, respectively.
SLID_2011_1S2 <- SLID_2011_1S1 %>% 
  filter(age_cat!=0 & educ_at!=99 & ml01v28!=30 & ml01v28!=99)
SLID_2011_2S2 <- SLID_2011_2S1 %>% 
  filter(age_cat!=0 & educ_at!=99 & ml11v28!=30 & ml11v28!=99)
# This line of code simply demonstrates that the two dataframes align perfectly once 
# we eliminate all the observations with the characteristics just described.
MATCH <- inner_join(SLID_2011_1S2,SLID_2011_2S2,
                    by=c("ecsex99","age_cat","educ_at","ml01v28"="ml11v28"))

# Now we prepare a series of tables that calculate the number and proportion of 
# individuals who began the year in the employed and unemployed states who then ended
# the year in one of those states. Each table examines how these values vary according
# to the combination of individual characteristics (i.e., gender, age, and education)
# we consider. The last of these is the one we ultimately use in the analysis.
SLID_2011_1S3_TOT <- SLID_2011_1S2 %>% group_by(ml01v28) %>% 
  summarise(Employed_12=sum(`Employed (DEC)`),Unemployed_12=sum(`Unemployed (DEC)`))
SLID_2011_1S3_TOT <- SLID_2011_1S3_TOT %>% 
  mutate(Emp_12Pr=Employed_12/(Employed_12+Unemployed_12)) %>%
  mutate(Unemp_12Pr=Unemployed_12/(Employed_12+Unemployed_12))
SLID_2011_1S3_Gen <- SLID_2011_1S2 %>% group_by(ecsex99,ml01v28) %>% 
  summarise(Employed_12=sum(`Employed (DEC)`),Unemployed_12=sum(`Unemployed (DEC)`))
SLID_2011_1S3_Gen <- SLID_2011_1S3_Gen %>% 
  mutate(Emp_12Pr=Employed_12/(Employed_12+Unemployed_12)) %>%
  mutate(Unemp_12Pr=Unemployed_12/(Employed_12+Unemployed_12))
SLID_2011_1S3_Edu <- SLID_2011_1S2 %>% group_by(educ_at,ml01v28) %>% 
  summarise(Employed_12=sum(`Employed (DEC)`),Unemployed_12=sum(`Unemployed (DEC)`))
SLID_2011_1S3_Edu <- SLID_2011_1S3_Edu %>% 
  mutate(Emp_12Pr=Employed_12/(Employed_12+Unemployed_12)) %>%
  mutate(Unemp_12Pr=Unemployed_12/(Employed_12+Unemployed_12))
SLID_2011_1S3_Age <- SLID_2011_1S2 %>% group_by(age_cat,ml01v28) %>% 
  summarise(Employed_12=sum(`Employed (DEC)`),Unemployed_12=sum(`Unemployed (DEC)`))
SLID_2011_1S3_Age <- SLID_2011_1S3_Age %>% 
  mutate(Emp_12Pr=Employed_12/(Employed_12+Unemployed_12)) %>%
  mutate(Unemp_12Pr=Unemployed_12/(Employed_12+Unemployed_12))
SLID_2011_1S3_Gen_Edu <- SLID_2011_1S2 %>% group_by(ecsex99,educ_at,ml01v28) %>% 
  summarise(Employed_12=sum(`Employed (DEC)`),Unemployed_12=sum(`Unemployed (DEC)`))
SLID_2011_1S3_Gen_Edu <- SLID_2011_1S3_Gen_Edu %>% 
  mutate(Emp_12Pr=Employed_12/(Employed_12+Unemployed_12)) %>%
  mutate(Unemp_12Pr=Unemployed_12/(Employed_12+Unemployed_12))
SLID_2011_1S3_Gen_Age <- SLID_2011_1S2 %>% group_by(ecsex99,age_cat,ml01v28) %>% 
  summarise(Employed_12=sum(`Employed (DEC)`),Unemployed_12=sum(`Unemployed (DEC)`))
SLID_2011_1S3_Gen_Age <- SLID_2011_1S3_Gen_Age %>% 
  mutate(Emp_12Pr=Employed_12/(Employed_12+Unemployed_12)) %>%
  mutate(Unemp_12Pr=Unemployed_12/(Employed_12+Unemployed_12))
SLID_2011_1S3_Edu_Age <- SLID_2011_1S2 %>% group_by(educ_at,age_cat,ml01v28) %>% 
  summarise(Employed_12=sum(`Employed (DEC)`),Unemployed_12=sum(`Unemployed (DEC)`))
SLID_2011_1S3_Edu_Age <- SLID_2011_1S3_Edu_Age %>% 
  mutate(Emp_12Pr=Employed_12/(Employed_12+Unemployed_12)) %>%
  mutate(Unemp_12Pr=Unemployed_12/(Employed_12+Unemployed_12))
SLID_2011_1S3_ALL <- SLID_2011_1S2 %>% group_by(ecsex99,age_cat,educ_at,ml01v28) %>% 
  summarise(Employed_12=sum(`Employed (DEC)`),Unemployed_12=sum(`Unemployed (DEC)`))
SLID_2011_1S3_ALL <- SLID_2011_1S3_ALL %>% 
  mutate(Emp_12Pr=Employed_12/(Employed_12+Unemployed_12)) %>%
  mutate(Unemp_12Pr=Unemployed_12/(Employed_12+Unemployed_12))
SLID_2011_2S3_ALL <- SLID_2011_2S2 %>% group_by(ecsex99,age_cat,educ_at,ml11v28) %>% 
  summarise(Employed_12=sum(`Employed (DEC)`),Unemployed_12=sum(`Unemployed (DEC)`))
SLID_2011_2S3_ALL <- SLID_2011_2S3_ALL %>% 
  mutate(Emp_12Pr=Employed_12/(Employed_12+Unemployed_12)) %>%
  mutate(Unemp_12Pr=Unemployed_12/(Employed_12+Unemployed_12))

# Next, we convert the values in the last two columns of each dataframe into 16x2 matrices
# (we cannot carry out matrix multiplication using dataframes).
SLID_2011_1S4 <- as.matrix(SLID_2011_1S3_ALL[,7:8])
SLID_2011_2S4 <- as.matrix(SLID_2011_2S3_ALL[,7:8])

# This for loop conducts a sequence of matrix multiplications to predict the likelihood
# of being in each employment state in January 2012, by assuming state transition
# probabilities in December 2011 are similar to those in November 2011.
SLID_2011_S5 <- matrix(0,nrow=16,ncol = 2)
for (i in seq(1,15,2)){
  SLID_2011_S5[i:(i+1),] <- SLID_2011_1S4[i:(i+1),1:2] %*% SLID_2011_2S4[i:(i+1),1:2]
}

# In this step we transform the calculated figures into a dataframe, merge them with
# the data frames prepared in Step 3, and rename the columns so they are more easily
# interpreted.
SLID_2011_S5 <- as.data.frame(SLID_2011_S5)
SLID_2011_S6 <- bind_cols(SLID_2011_1S3_ALL,SLID_2011_2S3_ALL[,5:8],SLID_2011_S5)
SLID_2011_S6 <- SLID_2011_S6 %>% 
  rename(Emp_01.12=Employed_12,Unemp_01.12=Unemployed_12,EmpPr_01.12=Emp_12Pr,
         UnempPr_01.12=Unemp_12Pr,Emp_11.12=Employed_121,Unemp_11.12=Unemployed_121,
         EmpPr_11.12=Emp_12Pr1,UnempPr_11.12=Unemp_12Pr1,EmpPr_01.01=V1,
         UnempPr_01.01=V2)

# Save results ------------------------------------------------------------

save(SLID_2011_S6, file = paste(RDA, Filename, LM, ".Rda", sep = ""))