# R Code for Analysis of 2011 National Household Survey PUMF
# Created: 09/01/18; Last Modified: 04/20/19

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

# NOTE 042019: Changes in the RSTUDIOAPI package appear to have affected the
# getSourceEditorContext()$path operation, causing issues when referencing
# project components. Ultimately, I determined that by setting the default
# working directory to the DIABETES_MODEL subfolder, it is possible to
# eliminate a lot of unnecessary complexity, since the path to all project
# components can now defined relative to the output generated by the 
# HERE function. I will update the User Guide to communicate the need for
# this added step to prospective users. 
# 
# 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 the script to use tax rates for 2018.
# 
# NOTE 060718: Migrated much of the content in this script, including derivation of
# retirement income and Pharmacare deductibles/out-of-pocket pharmaceutical expenditures,
# to RWD&TP_Matrices_Construction_MMDDYY.
# 
# NOTE 011918: Code previously used wages and salaries as the basis for calculations.
# This new code uses employment income, which appears to also include income from
# self-employment. It should be noted that we do not account for the fact that
# self-employed individuals are required to pay both the employer and employee
# component of the CPP contribution.
# 
# NOTE 051518: Code revised to draw from 2011 NHS PUMF, rather than summary CANSIM
# tables. The rationale for this is twofold: 
#   1) No pre-processing of the data is required; users accessing the PUMF can
#   immediately carry out the calculations in the code.
#   2) Using the PUMF enables calculation of income quintiles, which is utilized in the
#   DPoRT 2.0 risk equations.
# 
# Analogous to the 2011 NHS data, OAS & GIS benefits are calculated from a dataset
# published by the Government of Canada at https://open.canada.ca/en, as opposed to
# the PDF document used in the previous version of this code. This, again, promotes
# transparency and reproduceability, and ensures the calculations can be implemented
# without any need for pre-processing on the user's behalf.

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

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

# Call requisite libraries
library(dplyr); library(tidyr); library(reshape2); library(tibble); library(rvest)
library(pdftools); library(magrittr); library(jsonlite); library(httr);
library(data.table); library(lubridate); library(here); library(openxlsx)

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

# Preparing the Data ------------------------------------------------------

# First, we obtain CPI data for the period from 2000 to 2017 using the new
# Statistics Canada API. I could not find a straightforward way to extract information
# directly from the website; however, I WAS able to download the desired tables in CSV
# format using the API, unzip the folder in which it is contained, and then manipulate
# the dataset to get the information I wanted.
# https://www.r-bloggers.com/accessing-apis-from-r-and-a-little-r-programming/
# https://www150.statcan.gc.ca/t1/wds/

# 090618: "18100005" is a reference to "Consumer Price Index, annual average, not 
# seasonally adjusted" - Table: 18-10-0005-01 (formerly: CANSIM 326-0021). This had not
# been updated for 2018 as of the date this note was inserted.
RAW_DATA <- GET("https://www150.statcan.gc.ca/t1/wds/rest/getFullTableDownloadCSV/18100005/en")
FILE_LOC <- rawToChar(RAW_DATA$content)
STRING.00 <- nchar(FILE_LOC)
STRING.01 <- unlist(gregexpr('https', FILE_LOC))
STRING.02 <- substr(FILE_LOC, STRING.01, STRING.00-2)
# https://stackoverflow.com/questions/3053833/using-r-to-download-zipped-data-file-extract-and-import-data
temp <- tempfile()
download.file(STRING.02, temp)
CPI <- read.csv(unz(temp, "18100005.csv"))
unlink(temp)

CPI <- CPI %>%
  rename("Year" = "ï..REF_DATE", "PRODUCTS" = "Products.and.product.groups", 
         "CPI_2002" = "VALUE") %>%
  filter(Year %in% 2002:2017, GEO == "Canada", UOM_ID == 17, 
         PRODUCTS == "All-items") %>%
  select(Year, CPI_2002)

# Now we rebase the data to 2010 (the year from which the most recent NHS data is
# available).
Base_2010 <- CPI %>% filter(Year==2010) %>% select(CPI_2002)
Base_2010 <- data.matrix(Base_2010)[1,1]
CPI <- CPI %>% mutate(CPI_2010 = CPI_2002/Base_2010)
# Rebasing the CPI enables us to derive a deflator that can be used to convert the NHS
# data into 2017 dollars.
Deflator_2017 <- CPI %>% filter(Year==2017) %>% select(CPI_2010)
Deflator_2017 <- data.matrix(Deflator_2017)[1,1]

# The other data source is the 2011 NHS PUMF (99-014-X2011041). Note that the user will
# need to obtain this data source themselves (i.e., it will not be uploaded to MSpace,
# at least in part because this is not necessary).
# http://www12.statcan.gc.ca/nhs-enm/2011/dp-pd/prof/index.cfm?Lang=E
# http://www.dummies.com/programming/r/how-to-use-read-csv-to-import-data-in-r/
Manitoba_NHS_v1 <- 
  fread(paste0(Data_Path, "NHS-99M001X-E-2011-pumf-individuals_F1.csv", sep = ""))

# Extract Manitoba records
Manitoba_NHS_v1 <- Manitoba_NHS_v1 %>% 
  filter(PR == 46, WRKACT == 12, AGEGRP %in% c(9:21), HDGREE %in% c(1:13), 
         EMPIN != 8888888, EMPIN > 1) %>%
  select(PR, WRKACT, AGEGRP, SEX, HDGREE, EMPIN) %>%
  mutate(AGEGRP.01 = case_when(
    AGEGRP %in% 9:10 ~ "25_34",
    AGEGRP %in% 11:12 ~ "35_44",
    AGEGRP %in% 13:14 ~ "45_54",
    AGEGRP %in% 15:16 ~ "55_64",
    AGEGRP %in% 17:21 ~ "65_")) %>%
  mutate(HDGREE.01 = case_when(
    HDGREE %in% 1 ~ "1. No certificate, diploma or degree",
    HDGREE %in% 2 ~ "2. High school diploma or equivalent",
    HDGREE %in% 3:8 ~ "3. Postsecondary certificate or diploma below bachelor level",
    HDGREE > 8 ~ "4. University certificate, diploma or degree at bachelor level or above")) %>%
  mutate(Sex = case_when(SEX == 1 ~ "Female", SEX == 2 ~ "Male")) %>%
  mutate(EMPIN.01 = EMPIN*Deflator_2017)

Manitoba_NHS_QUIN <- Manitoba_NHS_v1 %>%
  group_by(Sex) %>% summarise(TOPQUIN_EMPIN = quantile(EMPIN.01, probs = 0.80))

Manitoba_NHS_v2 <- Manitoba_NHS_v1 %>%
  group_by(AGEGRP.01, Sex, HDGREE.01) %>%
  summarise(MED_EMPIN = median(EMPIN.01)) %>%
  left_join(., Manitoba_NHS_QUIN, by = "Sex") %>%
  mutate(IsTopQuin = MED_EMPIN > TOPQUIN_EMPIN)

# For the analysis of employment income, we take a subset of the data contained in 
# Manitoba_NHS_v2 and manipulate some columns (i.e., so that we can ultimately prepare
# separate tabulations for individuals who are full- and part-time).
Emp_Income <- Manitoba_NHS_v2 %>% 
  rename("Age_Cat" = "AGEGRP.01", "Highest_Degree" = "HDGREE.01", "Med_Emp_Inc_2017" = 
     "MED_EMPIN") %>%
  select(Age_Cat, Sex, Highest_Degree, Med_Emp_Inc_2017, IsTopQuin) %>%
  rename(Med_Emp_Inc_2017.FT = Med_Emp_Inc_2017) %>%
  mutate(Med_Emp_Inc_2017.PT = Med_Emp_Inc_2017.FT/2)

RPT_Emp_Income <- 
  dcast(Emp_Income[, c(1:4)], Highest_Degree + Sex ~ Age_Cat, 
        value.var = "Med_Emp_Inc_2017.FT") %>%
  arrange(., Highest_Degree, desc(Sex))

wb = createWorkbook(); sht = addWorksheet(wb, "Gross_Income")
dat = RPT_Emp_Income; writeData(wb, sht, dat, colNames = T)
saveWorkbook(wb, paste0(Model_Path, "/Report_Tables/", "Gross_Income",
                        LM, ".xlsx"), overwrite = T)

save(Emp_Income, file = paste0(RDA, Filename, LM, "_1", ".Rda"))

# Calculating Employment & Retirement Income --------------------------------------------------------

# Tax rates reflect combined federal and provincial taxes. The pertinent values for
# 2018 are drawn from: https://www.ey.com/Publication/vwLUAssets/Tax-Rates-Manitoba-2018/%24FILE/Tax-Rates-Manitoba-2018.pdf
Tax_Rates <- pdf_text(paste("https://www.ey.com/Publication/vwLUAssets/",
                            "Tax-Rates-Manitoba-2018/%24FILE/",
                            "Tax-Rates-Manitoba-2018.pdf", sep = ""))
Tax_Rates <- Tax_Rates[1]

# https://www.regular-expressions.info/
# https://stackoverflow.com/questions/10617702/remove-part-of-string-after
# https://stackoverflow.com/questions/15895050/using-gsub-to-extract-character-string-before-white-space-in-r
Tax_Rates %<>%
  # gsub("\r\n", "", .) %>%
  # Eliminate all alphabetical characters
  gsub('([[:alpha:]])', "",.) %>%
  # Eliminate percentage signs
  gsub("%", "",.) %>%
  # Convert all spaces consisting of two or more consecutive blank characters to a
  # single blank character
  gsub(" {2,}", " ",.) %>%
  # Eliminate everything from the string "\r\n1" up until the end of the string.
  gsub("\r\n1.*$", "",.) %>%
  # Eliminate everything from the beginning of the string (denoted by "^"), up to and
  # including "4\r\n".
  gsub("^.*4\r\n", "",.) %>%
  # Eliminate the first two characters of the string.
  sub(".{2}", "",.) %>%
  # Replace commas with periods (for some reason the original document uses commas in
  # the place of a decimal point).
  gsub(",", "",.) %>%
  # Replace all currency signs with a space. The "\\" is necessary to ensure the 
  # character is actually interpreted as a dollar sign, and not as the end of the
  # string (see above for an example of how the dollar sign is applied to serve this
  # alternative function).
  gsub(" \\$ ", " ",.)

# The substring "\r\n" now separates each row of the original table, so we'll split the
# whole string using this as a delimiter.
Tax_Rates <- unlist(strsplit(Tax_Rates, "\r\n"))
# The last row contains a blank where there ought to be (or where it would be 
# convenient for there to be) a dash, so we insert one now.
Tax_Rates[9] <- sub("3 7", "3 - 7", Tax_Rates[9])
# Eliminate the first character in every string.
Tax_Rates <- as.character(sapply(Tax_Rates, FUN = function(x) gsub("^.{1}", "",x)))

# The spaces between the numbers in the first three columns is problematic, because R
# has no way of knowing these demarcate parts of the SAME NUMBERS (in most documents,
# this would be done with a comma). Thus, if we're going to use spaces as a deliminter
# (i.e., to distinguish between numbers and not between their parts), we'll need to
# address this.
for (i in 1:9){
  # Unfortunately, the values in each row of the first three columns of the original
  # table typically contain different numbers of characters, so we need to address
  # these on a case by case basis.
  Tax_Rates[i] <- case_when(
    # NOTE: "\d" refers to a digit, and the number in {} denotes how many characters we
    # are looking for (see: http://stat545.com/block022_regular-expression.html). It
    # appears that "-" is not interpreted as a digit. 
    # 
    # R interprets the expression following the first column as a series of references 
    # to distinct items presented before the comma. Hence, (\\d{1}) is Item #1 (i.e.,
    # "\\1"), (\\d{3}) is Item #2 (i.e., "\\2"), and so on. Thus, the expression
    # indicates that we want to concatenate the first and second items and the third
    # fourth items. The traditional (for lack of a better word) parentheses (i.e., "("
    # and ")") presumably exist to demarcate the numbered items.
    i == 1 ~ gsub('(\\d{1}) (\\d{3})','\\1\\2', Tax_Rates[i]),
    i %in% 2:3 ~ gsub('(\\d{1}) (\\d{3}) (\\d{2}) (\\d{3})','\\1\\2 \\3\\4', Tax_Rates[i]),
    i %in% 4:5 ~ gsub('(\\d{2}) (\\d{3}) (\\d{2}) (\\d{3}) (\\d{1}) (\\d{3})','\\1\\2 \\3\\4 \\5\\6', Tax_Rates[i]),
    i == 6 ~ gsub('(\\d{2}) (\\d{3}) (\\d{2}) (\\d{3}) (\\d{2}) (\\d{3})','\\1\\2 \\3\\4 \\5\\6', Tax_Rates[i]),
    i == 7 ~ gsub('(\\d{2}) (\\d{3}) (\\d{3}) (\\d{3}) (\\d{2}) (\\d{3})','\\1\\2 \\3\\4 \\5\\6', Tax_Rates[i]),
    i == 8 ~ gsub('(\\d{3}) (\\d{3}) (\\d{3}) (\\d{3}) (\\d{2}) (\\d{3})','\\1\\2 \\3\\4 \\5\\6', Tax_Rates[i]),
    i == 9 ~ gsub('(\\d{3}) (\\d{3}) - (\\d{2}) (\\d{3})','\\1\\2 - \\3\\4', Tax_Rates[i])
  )
}

# Now that each string is structured the way we want it, we need to split up the
# contents by column, where columns are delimited simply by the presence of a space.
Tax_Rates.02 <- lapply(Tax_Rates, FUN = function(x) unlist(strsplit(x, split = " ")))
# Generate column names.
Tax_Names <- c("Lower.Limit", "Upper.Limit", "Basic.Tax", "MR", "MR.Elg.Dividend",
               "MR.Other.Dividend", "MR.Capital.Gains")

Tax_Rates.03 <- as.data.frame(matrix(unlist(Tax_Rates.02), nrow = 9, byrow = TRUE))
names(Tax_Rates.03) <- Tax_Names

Tax_Rates.03 <- as.data.frame(apply(Tax_Rates.03, 2, FUN = function(x) as.numeric(as.character(x))))
Tax_Rates.03[,4:7] <- Tax_Rates.03[,4:7]/100
Tax_Rates.03$Upper.Limit[9] <- 1000000
Tax_Rates.03[is.na(Tax_Rates.03)] <- 0

# To simplify the calculations, we calculate the amount of taxes an individual would
# pay if they had incomes at the upper range of each tax bracket (i.e., the marginal
# rate, multiplied by the size of the bracket). Then, we calculate cumulative taxes
# paid.
Tax_Rates <- Tax_Rates.03 %>% 
  mutate(Tax_Amt = MR * (Upper.Limit-Lower.Limit)) %>%
  mutate(Cum_Tax_Amt = lag(cumsum(Tax_Amt)))

# When we apply the lag, the first value is null (i.e., N/A). We address this by
# replacing the value in this cell by 0.
Tax_Rates$Cum_Tax_Amt[is.na(Tax_Rates$Cum_Tax_Amt)==TRUE] <- 0

save(Tax_Rates, file = paste0(RDA, Filename, LM, "_2", ".Rda"))
wb = createWorkbook(); sht = addWorksheet(wb, "Tax_Rates")
dat = Tax_Rates; writeData(wb, sht, dat, colNames = T)
saveWorkbook(wb, paste0(Model_Path, "/Report_Tables/", "Tax_Rates",
                        LM, ".xlsx"), overwrite = T)