Project-OLS.Rmd

---
title: "Estimation strategy"
author: "Sumit Meghlani 64801343"
date: "3/25/2022"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(tibble)
library(dplyr)
library(stargazer)
library(tidyverse)
library(fixest)
library(etable)
library(modelsummary)
library(kableExtra)
library(FactoMineR)
```

```{r functions}
#' @param df dataframe (decide if rows have been randomized)
#'    randomized outside of function if mix R/python methods
#' @param K number of folds for cross-validation
#' @param nperfmeas number of performance measures used to compare methods
#' @param seed random number of seed for reproducibility
#' @return matrix of performance measures by fold, vector of averages
crossValidationCont = function(df,K,nperfmeas=6,seed)
{ set.seed(seed)
  # apply transformations
  df <- transform(df)
  
  n = nrow(df)
  nhold = round(n/K) # size of holdout set 
  iperm = sample(n)
  comperf50meas = matrix(0,K,nperfmeas)
  comperf80meas = matrix(0,K,nperfmeas)
  resperf50meas = matrix(0,K,nperfmeas)
  resperf80meas = matrix(0,K,nperfmeas)
  
  for(k in 1:K)
  { indices = (((k-1)*nhold+1):(k*nhold))
    if( k==K ) indices = (((k-1)*nhold+1):n)
    indices = iperm[indices]
    # split between commercial and residential 
    commdf <- ClassSplit(df)[[1]]%>% medianImp()
    resdf <- ClassSplit(df)[[2]]%>% medianImp()
    
    # commercial
    commtrain = commdf[-indices,]
    commholdout = commdf[indices,]
    
    PredfestcomModel <- feols(site_eui~minavg+maxavg+energy_star_rating+floor_area+facility_type+Year_Factor+State_Factor+cooling_degree_days+heating_degree_days+precipitation_inches+days_above_110F+days_below_30F,commtrain)
    
    restrain = resdf[-indices,] 
    resholdout = resdf[indices,] 
    PredfestresModel <- feols(site_eui~minavg+maxavg+energy_star_rating+facility_type+Year_Factor+State_Factor+cooling_degree_days+heating_degree_days+precipitation_inches+days_above_110F+days_below_30F,restrain)
    #pred = myPredict(obj, newdata=holdout)  # could include 50% and 80% intervals
    
    comperf50meas[k,] = probCheck(PredfestcomModel,commholdout,0.51)$summary
    comperf80meas[k,] = probCheck(PredfestcomModel,commholdout,0.80)$summary
    resperf50meas[k,] = probCheck(PredfestresModel,resholdout,0.51)$summary
    resperf80meas[k,] = probCheck(PredfestresModel,resholdout,0.80)$summary
      #perfMeas(pred, holdout$y)
  }
  # 50
  comavgperf50meas = apply(comperf50meas,2,mean)
  resavgperf50meas = apply(resperf50meas,2,mean)
  # 80 
  comavgperf80meas = apply(comperf80meas,2,mean)
  resavgperf80meas = apply(resperf80meas,2,mean)
  
  list(comperfmeas50byfold=comperf50meas,resperfmeas50byfold = resperf50meas, comavgperf50meas=comavgperf50meas,resavgperf50meas = resavgperf50meas,comperfmeas80byfold=comperf80meas,resperfmeas80byfold = resperf80meas, comavgperf80meas=comavgperf80meas,resavgperf80meas = resavgperf80meas)
}


# prob matrix
probCheck <- function(model,df,level){
  problist <- predict(model,df,interval='predi',level=level)[,-2] %>% na.omit()
  return(intervalScore(problist,na.omit(df$site_eui),level))
}
```


```{r functions}
statebins <- function(df,ntrain){
  weather= rep("M",ntrain)
  weather[df$State_Factor == 'State_4']='C';weather[df$State_Factor == 'State_6']='C';weather[df$State_Factor == 'State_1']='H';weather[df$State_Factor == 'State_11']='H';
  weather = as.factor(weather); df$State_Factor = weather
  return(df)
}

transform <- function(df,subset = TRUE){
  if (subset==TRUE){
    df <- df
    ntrain <- nrow(df)/2
    df <- df[sample(nrow(df),ntrain),]
    df <-  statebins(df,ntrain)
    }
    # transformations
    df$floor_area <- log(df$floor_area)
    df$site_eui <- log(df$site_eui)
    
    # days_data
    daysDatNames <- df %>% select(contains('days')) %>% names()
    
    # creating auxillarydf
    #auxDf <- df %>% select(c(daysDatNames,'avg_temp',"cooling_degree_days","heating_degree_days","precipitation_inches","snowfall_inches","snowdepth_inches","building_class"))
    
    # removing missing values
    
    #df <- df %>% select(-c('id','direction_peak_wind_speed' ,'max_wind_speed' ,'days_with_fog',"cooling_degree_days","heating_degree_days")) 
    # removed percip,snow,days
    
    #df <- df %>% select(-names(df)[45:57])
    # min and max names
    minNames <- df %>% select(contains("min")) %>% names()
    maxNames <- df%>% select(contains("max")) %>% names()
    avgNames <- df%>% select(contains("_avg_")) %>% names()
    
    # create a minimum and maximum averagese.
    df <- df %>% mutate(minavg = (january_min_temp+february_min_temp+march_min_temp+october_min_temp+november_min_temp+december_min_temp)/4)
    df <- df %>% mutate(maxavg = (april_max_temp+may_max_temp+june_max_temp+july_max_temp+august_max_temp+september_max_temp)/4)
    
    #return(list(df,auxDf))
    return(df)
  
}

medianImp <- function(df){
  df$year_built[is.na(df$year_built)]= summary(df$year_built)[[3]]
  df$energy_star_rating[is.na(df$energy_star_rating)] = summary(df$energy_star_rating)[[3]]
  return(df)
}

ClassSplit <- function(df){
  commdf <- df %>% filter(building_class=='Commercial')%>% select(-'building_class')
  resDf <- df %>% filter(building_class=="Residential") %>% select(-'building_class')
  return(list(commdf,resDf))
}

plot_results <- function(df,model){
  preds <- predict(model,df)
  plot(density(preds));
  lines(density(df$site_eui),col=2,type='h',main='predictions vs actual');
}

intervalScore = function(predObj,actual,level)
{ n = nrow(predObj)
  alpha = 1-level
  ilow = (actual<predObj[,2])  # underestimation
  ihigh = (actual>predObj[,3]) # overestimation
  sumlength = sum(predObj[,3]-predObj[,2]) # sum of lengths of prediction intervals
  sumlow = sum(predObj[ilow,2]-actual[ilow])*2/alpha
  sumhigh = sum(actual[ihigh]-predObj[ihigh,3])*2/alpha
  avglength = sumlength/n
  IS = (sumlength+sumlow+sumhigh)/n # average length + average under/over penalties
  browser()
  cover = mean(actual>= predObj[,2] & actual<=predObj[,3])
  summ = c(level,avglength,IS,cover) 
  # summary with level, average length, interval score, coverage rate
  imiss = which(ilow | ihigh)
  list(summary=summ, imiss=imiss)
}
```

## loading in the dataset
```{r}
# load the data
bigdf <- read.csv('train.csv')

# use this to omit na values
# bigdf <- na.omit(bigdf)
ntrain = nrow(bigdf)
df = bigdf[sample(nrow(bigdf),ntrain),]
# to use the full dataset

# log transformation for floor_area
# going to use max of summer months [Apr - Sep], min of winter months [Oct - Mar]
```


```{r}
df <- transform(df=bigdf)
# ClassSplit the main df
commDf <- ClassSplit(df)[[1]] %>% medianImp()
resDf <- ClassSplit(df)[[2]] %>% medianImp()
```

### Main Model regression equation

**$e_{t,J} = \alpha_{t,J} + \beta_1 minavg_{t,J} + \beta_2 maxavg_{t,J} + \beta_3 ESA_{t,J} +\beta_4 yb_{t,J} +\beta_5 elev_{t,J}+\beta_6 floar_{t,J}+\beta_7 factyp_{t,J}+ \epsilon_{t,J}$**

# running fest
Fest is a package in R used for fixed effects estimation.

```{r}
#running fixest
#-- cant be used for predictions as standard errors harder to compute with fixed effects
festcomModel <- feols(site_eui~minavg+maxavg+energy_star_rating+facility_type|Year_Factor+State_Factor|avg_temp~cooling_degree_days+heating_degree_days+precipitation_inches+days_above_110F+days_below_30F,commDf)

festresModel <- feols(site_eui~+minavg+maxavg+energy_star_rating+facility_type|Year_Factor+State_Factor|avg_temp~cooling_degree_days+heating_degree_days+precipitation_inches+days_above_110F+days_below_30F,resDf)
```

```{r}
# - plotting density predictions
plot_results(commDf,festcomModel)
plot_results(resDf,festresModel)
```

# interval score implementation
```{r}
# running models to predict
PredfestcomModel <- feols(site_eui~minavg+maxavg+energy_star_rating+facility_type+Year_Factor+State_Factor|avg_temp~cooling_degree_days+heating_degree_days+precipitation_inches+days_above_110F+days_below_30F,commDf)

PredfestresModel <- feols(site_eui~+minavg+maxavg+energy_star_rating+facility_type+Year_Factor+State_Factor+cooling_degree_days+heating_degree_days+precipitation_inches+days_above_110F+days_below_30F,resDf)

```


```{r}
# 0.5 predictions
comModel50 <- probCheck(PredfestcomModel,commDf,0.51) 
resModel50 <- probCheck(PredfestresModel,resDf,0.51) 
#0.8 predictions
comModel80 <- probCheck(PredfestcomModel,commDf,0.95)
resModel80 <- probCheck(PredfestresModel,resDf,0.95)

REcomreg <- rbind(comModel50$summary,comModel80$summary)
colnames(REcomreg)=c("level","avgleng","IS","cover")
print(REcomreg)
REresreg <- rbind(resModel50$summary,resModel80$summary)
colnames(REresreg)=c("level","avgleng","IS","cover")
print(REresreg)
```





```{r residential feeffects echo=FALSE}
# running fixed effects one by one
feList = list()
all_FEs = c("State_Factor","Year_Factor")
for (i in 0:2){
  feList[[i+1]] <- feols(log(site_eui)~+minavg+maxavg+energy_star_rating+year_built+ELEVATION+floor_area+i(facility_type)+avg_temp,resDf,fixef = all_FEs[0:i])
}
etable(feList,title='Residential Regression Table')

cm <- c('(Intercept)'='Intercept','minavg'='Minimum Average (F)','maxavg'='Maximum Average(F)','energy_star_rating'='Energy Star Rating','year_built'='Year Built','floor_area'='Floor Area','avg_temp'='Average Temp (aux)')
tb <- modelsummary::modelsummary(feList,output="kableExtra",stars=TRUE,title = 'Regression Table',coef_map = cm,gof_omit = "BIC|Log.Lik.|R2 Adj.|R2 Within| R2 Pseudo|AIC")

```

#crossValidation
```{r}
crossValidationCont(bigdf,5,nperfmeas = 4,seed=123)
```

#Principal Component Analysis
Cannot do PRincipal Componenet Analysis for the life of me, if someone figures this out. Thanks in advance.
```{r}
#test <- tibble(commDf$facility_type)
#princomp(na.omit(commDf[,3]),na.action=na.exclude, cor = TRUE)
#testmfa = MFA(base=test,group=c(54),type = c("n"),ncp=5,name.group = c("facility_type"))
```

# see below the residential coefs plot

```{r residential coefficient estimate plots}
plotDf <- sort(resfestModel$coefficients)
qplot(x=plotDf,y = names(plotDf))+geom_bar(stat="identity", position="dodge")+labs(title='Residential Regression Coefficients',x = '%',y='Features')
```

```{r commercial fest for loop echo=FALSE}
feList = list()
all_FEs = c("State_Factor","Year_Factor")
for (i in 0:2){
  feList[[i+1]] <- feols(log(site_eui)~+minavg+maxavg+energy_star_rating+year_built+ELEVATION+floor_area+i(facility_type)+avg_temp,commDf,fixef = all_FEs[0:i])
}
etable(feList,drop='=',tex=TRUE,title='Commercial Regression Table')
```
Models to run interval score on
- festModel
- resfestModel
```{r}

plot(density(predict(festModel,commDf)))
lines(density(log(commDf$site_eui)),col=2,type='h')
```