---
title: "Data supporting \"Evaluating the Opportunities and Challenges of Domestic PV Installation in Saudi Arabia Based on Field Deployment in Jeddah\""
author: "Luke Blunden [lsb1@soton.ac.uk](mailto:lsb1@soton.ac.uk)"
date: "May 2025"
output:
  html_document:
    toc: true
    number_sections: true
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = FALSE)
knitr::opts_chunk$set(include = FALSE)
```

# Introduction
This R Markdown script is for reproducing the analysis and graphs from the journal paper:

> Alghamdi, A.; Blunden, L.S.; Alam, M.; Bahaj, A.S.; James, P.A.B. "Evaluating the Opportunities and Challenges of Domestic PV Installation in Saudi Arabia Based on Field Deployment in Jeddah". Energies 2025, 18, x. https://doi.org/10.3390/xxxxx

Included along with this file are the data files and a compiled html version of this script.  

# Load libraries

```{r the_libs, echo=TRUE, include=TRUE, warning=FALSE, error=FALSE, message=FALSE, results= 'hide'}
library(readxl)
library(dplyr)
library(magrittr)
library(ggplot2)
library(purrr)
library(lubridate)
library(tidyr)
library(plotly)
library(paletteer)
library(stargazer)
library(htmltools)
library(scales)
library(FinancialMath)
library(sjPlot)
library(kableExtra)
```

# Assign parameters

```{r params, include=TRUE, echo=TRUE}
# change filepath as appropriate
path_to_data <- ""
top_string <- as_labeller(c(`1` = "Villa 1", `2` = "Villa 2"))
TFA_1st <- 300 # first floor area, sq. m.
PV_area <- 70.6 # sq.m PV array
base_temp <- 18.0 # degC, for calculating CDDs
```

# Load data
## PV performance data

```{r get_PV_perf_data}
fname <- "JeddahPV.xlsx"
the_path <- paste0(path_to_data,fname)
the_cols <- read_xlsx(the_path,n_max=2,col_names = F,sheet="PVcombined")
for (i in 1:dim(the_cols)[1]) {
  for (j in 1:dim(the_cols)[2]){
    the_cols[1,j] %<>% gsub("([^|]*) \\|.*","\\1",.) %>% gsub("[-.: ]","_",.) 
  }
}
the_units <- unlist(the_cols[2,]) # unlist() extracts vector from object 
the_cols <- unlist(the_cols[1,])
the_data <- read_xlsx(the_path,sheet="PVcombined",skip = 2,col_names = the_cols,na = "n/a")
the_data$ptime <- strptime(the_data$Date_and_time,"%d.%m.%Y",tz = "Asia/Riyadh")
the_data$Date_and_time <- NULL # drop, don't need
the_data$Villa %<>% as.factor()
the_data %<>% mutate(
                      efrac=Consumed_directly/Consumption,
                      ifrac=Consumed_directly/(Insolation_energy*PV_area*(1/3600))
                      )
dat2019 <- read_xlsx(the_path,col_names = T,skip=1,sheet="dat2019")
dat2019_long <- pivot_longer(dat2019,Jan:Dec,names_to = "Month")
dat2019_long %<>% mutate(Month.ct=factor(Month, month.abb, ordered = TRUE))
```

## Billing data

```{r get_bill_data}
the_path_bill <- paste0(path_to_data,fname)
the_bill_data <- read_xlsx(the_path_bill,sheet="bill_data")
the_bill_data$Villa %<>% as.factor()
the_bill_data %<>% group_by(Villa) %>% mutate(cumConsKWh=cumsum(ConsKWh))
the_bill_data$BillMonth %<>% force_tz(tzone = "Asia/Riyadh")

kWhZeroVilla1 <- 
  the_bill_data$cumConsKWh[(the_bill_data$BillMonth == strptime("01.05.2022","%d.%m.%Y",tz = "Asia/Riyadh")) & (the_bill_data$Villa == 1)]
kWhZeroVilla2 <- 
  the_bill_data$cumConsKWh[(the_bill_data$BillMonth == strptime("01.05.2022","%d.%m.%Y",tz = "Asia/Riyadh")) & (the_bill_data$Villa == 2)]

the_bill_data$BillMonth %<>% +days(14) # approx centre in month for plotting
```

## Weather data from KAU 

```{r cp_met}
# car port met data
fname_cpmet <- "met.csv"
path_cpmet <- paste0(path_to_data,fname_cpmet)
cpmet <- read.csv(path_cpmet,stringsAsFactors = F)
cpmet$ptime <- strptime(cpmet$mt_time,"%Y-%m-%d %H:%M:%S",tz = "Asia/Riyadh")
cpmet %<>% arrange(desc(row_number()))
cpmet %<>% filter(ptime > as.POSIXct("2020-12-31 23:59:59", tz = "Asia/Riyadh"))

cpmet_daily_stats <- cpmet %>%
  group_by(ymd=format(ptime, format="%Y%m%d")) %>%
  summarise(
            maxta=max(temp_degC),
            minta=min(temp_degC),
            avta=mean(temp_degC),
            avta2=(maxta+minta)/2,
            medta=median(temp_degC),
            cdd=sum(if_else(temp_degC>18.0,(temp_degC-18.0)*(1/720),0),na.rm = TRUE),
            maxRH=max(RH_pc),
            minRH=min(RH_pc),
            avRH=mean(RH_pc),
            medRH=median(RH_pc),
            maxIr=max(Irr_Wpsqm),
            sumIr=PV_area*(1/12)*(1/1000)*(sum(Irr_Wpsqm,na.rm = TRUE)), # 5 min interval
            ptime=mean(ptime),
            nobs=n(),
            .groups = "drop",
            na.rm = TRUE
            )

cpmet_daily_stats %<>% mutate(cdd2=if_else(avta>18.0,(avta-18.0),0))

cpmet_daily_stats_clean <- cpmet_daily_stats %>%  mutate(nobs > 720*0.9)  # only keep days where 90% or more of data available
```

## Weather data from rooftop mounted weather station

```{r get_met_data}
met_path <- paste0(path_to_data,fname)
met_cols <- read_xlsx(met_path,n_max=2,col_names = F,sheet="met_villa1_5min")
for (i in 1:dim(met_cols)[1]) {
  for (j in 1:dim(met_cols)[2]){
    met_cols[1,j] %<>% gsub("([^|]*) \\|.*","\\1",.) %>% gsub("[-.: ]","_",.) 
  }
}
met_units <- unlist(met_cols[2,]) # unlist() extracts vector from object 
met_cols <- unlist(met_cols[1,])

# Villa 1
met_data1 <- read_xlsx(met_path,sheet="met_villa1_5min",skip = 2,col_names = met_cols,na = "n/a")
met_data1$ptime <- strptime(met_data1$Date_and_time,"%d.%m.%Y %H:%M",tz = "Asia/Riyadh")
met_data1$Date_and_time <- NULL # drop, don't need
dt1 <- int_diff(met_data1$ptime) # check for regular time series

# Villa 2
met_data2 <- read_xlsx(met_path,sheet="met_villa2_5min",skip = 2,col_names = met_cols,na = "n/a")
met_data2$ptime <- strptime(met_data2$Date_and_time,"%d.%m.%Y %H:%M",tz = "Asia/Riyadh")
met_data2$Date_and_time <- NULL # drop, don't need
dt2 <- int_diff(met_data2$ptime) # check for regular time series
```

# Aggregate PV data to monthly to compare with bill data

```{r aggPV}
monthly_sums <- the_data %>%
  group_by(Villa=Villa,ym=format(ptime, format="%Y%m")) %>%
  summarise(
            across(.cols = c(Energy:Energy_MPP2),.fns = ~sum(.x)/1000),
            across(.cols = c(Consumed_directly:Energy_to_grid),.fns = ~sum(.x)/1000),
            across(.cols = c(Insolation_energy),.fns = ~70.6*sum(.x)/1000/3600),
            across(.cols = c(ptime),.fns = mean),
            .groups = "drop"
            )

monthly_sumsV1 <- monthly_sums %>% filter(Villa==1) %>% mutate(cumEnergy_from_grid=cumsum(Energy_from_grid))
monthly_sumsV2 <- monthly_sums %>% filter(Villa==2) %>% mutate(cumEnergy_from_grid=cumsum(Energy_from_grid))


monthly_sums <- bind_rows(monthly_sumsV1,monthly_sumsV2)

monthly_sums %<>% mutate(cumEnergy_from_grid=if_else(Villa==1,cumEnergy_from_grid+kWhZeroVilla1,cumEnergy_from_grid+kWhZeroVilla2))

monthly_sums_long <- pivot_longer(monthly_sums,cols = c(Consumption,Energy_from_grid,Consumed_directly))

```

# Join the bill data to the aggregated electrical data

```{r joined_up}
bill_data_copy=data.frame(
                          ptime=the_bill_data$BillMonth,
                          ym=format(the_bill_data$BillMonth, format="%Y%m"),
                          ConsKWh=the_bill_data$ConsKWh,
                          Villa=the_bill_data$Villa
                          )

joined_up_data <- bind_rows(bill_data_copy,monthly_sums)
joined_up_data %<>% mutate(cons2=if_else(is.na(Consumption),ConsKWh/TFA_1st,Consumption/TFA_1st)) # per floor area

joined_up_data_long <- joined_up_data %>%
  pivot_longer(c(ConsKWh,Consumed_directly,))


```

# Aggregate weather data

```{r agg_met_day}
met_daily_stats1 <- met_data1 %>%
  group_by(ymd=format(ptime, format="%Y%m%d")) %>%
  summarise(
            maxta=max(Ambient_temperature),
            minta=min(Ambient_temperature),
            avta=mean(Ambient_temperature),
            medta=median(Ambient_temperature),
            cdd=sum(if_else(Ambient_temperature>18.0,(Ambient_temperature-base_temp)*(1/288),0),na.rm = TRUE),
            maxtpv=max(Module_temperature),
            mintpv=min(Module_temperature),
            avtpv=mean(Module_temperature),
            medtpv=median(Module_temperature),
            maxIr=max(Irradiation),
            sumIr=PV_area*(1/12)*(1/1000)*(sum(Irradiation,na.rm = TRUE)), # 5 min interval
            ptime=mean(ptime),
            .groups = "drop",
            na.rm = TRUE
            )

met_daily_stats1 %<>% mutate(cdd2=if_else(avta>base_temp,(avta-base_temp),0))

met_daily_stats2 <- met_data2 %>%
  group_by(ymd=format(ptime, format="%Y%m%d")) %>%
  summarise(
            maxta=max(Ambient_temperature),
            minta=min(Ambient_temperature),
            avta=mean(Ambient_temperature),
            medta=median(Ambient_temperature),
            cdd=sum(if_else(Ambient_temperature>base_temp,(Ambient_temperature-base_temp)*(1/288),0),na.rm = TRUE),
            maxtpv=max(Module_temperature),
            mintpv=min(Module_temperature),
            avtpv=mean(Module_temperature),
            medtpv=median(Module_temperature),
            maxIr=max(Irradiation),
            sumIr=PV_area*(1/12)*(1/1000)*(sum(Irradiation,na.rm = TRUE)), # 5 min interval
            ptime=mean(ptime),
            .groups = "drop",
            na.rm = TRUE
            )

met_daily_stats2 %<>% mutate(cdd2=if_else(avta>base_temp,(avta-base_temp),0))
```

# Composite temperature dataset

```{r met_monthly}
met_daily_combined <- bind_cols(select(met_daily_stats1,c(ymd:medtpv,cdd2)),select(met_daily_stats2,c(maxIr:ptime)))

compare_met <- inner_join(cpmet_daily_stats_clean,met_daily_combined,by="ymd",suffix = c(".cp",".fr"))

fr_start <- min(met_daily_combined$ptime)

#plot(compare_met$avta.cp,compare_met$avta.fr)
mod_cp_fr_ta_0 <- lm(avta.cp~avta.fr+0,data = compare_met)
mod_cp_fr_ta <- lm(avta.cp~avta.fr,data = compare_met)
summary(mod_cp_fr_ta)
confint(mod_cp_fr_ta)
BIC(mod_cp_fr_ta_0) # BIC lower with zero intercept
BIC(mod_cp_fr_ta) 
tab2 <- stargazer(mod_cp_fr_ta_0,ci=T,type = "html")

calib_ta <- mod_cp_fr_ta_0$coefficients["avta.fr"] # calibration coefficient

cpmet_daily_stats_clean %<>% mutate(avta.calib = avta/calib_ta)

ta_composite <- bind_rows(select(
                                  filter(cpmet_daily_stats_clean,ptime<fr_start),
                                  c(ptime,avta=avta.calib)
                                  ),
                          select(met_daily_combined,c(ptime,avta))
                          ) 

ta_composite %<>% mutate(after=ptime>=fr_start)
ta_composite %<>% mutate(cdd=if_else(avta>base_temp,(avta-base_temp),0))


met_monthly_stats <- met_daily_combined %>% filter(ptime<"2024-06-03") %>%
  group_by(ym=format(ptime - days(2), format="%Y%m")) %>%
  summarise(
            cdd=sum(cdd),
            sumIr=sum(sumIr),
            ptime=mean(ptime),
            nobs=n(),
            .groups = "drop"
            )

Composite_met_monthly_stats <- ta_composite %>% filter(ptime<"2024-06-01") %>%
  group_by(ym=format(ptime, format="%Y%m")) %>%
  summarise(
            cdd=sum(cdd,na.rm = T),
            avta=mean(avta,na.rm = T),
            ptime=mean(ptime),
            nobs=n(),
            .groups = "drop"
            )

Composite_met_monthly_stats %<>% mutate(cdd.adj = cdd*days_in_month(ptime)/nobs)

bill_wide <- pivot_wider(bill_data_copy,names_from = Villa,names_prefix = "Villa_",names_sep = ".",values_from = ConsKWh)

composite_met_bill <- inner_join(Composite_met_monthly_stats,joined_up_data,by = c("ym"))
composite_met_bill %<>% mutate(after=ptime.x > "2022-05-31")

mod_cmbV1_l <- composite_met_bill %>% filter(Villa==1) %>% lm(cons2~cdd.adj, .)
sjPlot::tab_model(mod_cmbV1_l, vcov.fun = "HC3", show.se = TRUE)
summary(mod_cmbV1_l)
mod_cmbV2_l <- composite_met_bill %>% filter(Villa==2) %>% lm(cons2~cdd.adj, .)
summary(mod_cmbV2_l)
mod_cmbV1_q <- composite_met_bill %>% filter(Villa==1) %>% lm(cons2~cdd.adj+I(cdd.adj^2), .)
summary(mod_cmbV1_q)
mod_cmbV2_q <- composite_met_bill %>% filter(Villa==2) %>% lm(cons2~cdd.adj+I(cdd.adj^2), .)
summary(mod_cmbV2_q)
mod_cmbV2_nl <- lm(cons2~cdd.adj*after,data=composite_met_bill)

composite_met_bill %<>% mutate(
  mod_cons_l=case_when(
    Villa == 1 ~ predict(mod_cmbV1_l,.),
    Villa == 2 ~ predict(mod_cmbV2_l,.),
    )
  )

met_vs_cons <- inner_join(met_monthly_stats,monthly_sumsV1,by = c("ym"))
met_vs_cons %<>% inner_join(y=monthly_sumsV2,by = c("ym"),suffix = c("V1","V2"))
#met_vs_cons_long <- met_vs_cons %>% pivot_longer(c(ConsKWh,Consumed_directly))

mod_consV1cdd <- lm(ConsumptionV1~cdd,data = met_vs_cons)
mod_consV2cdd <- lm(ConsumptionV2~cdd,data = met_vs_cons)
summary(mod_consV1cdd)
confint(mod_consV1cdd)
tab1 <- stargazer(mod_cmbV1_l,mod_cmbV2_l,ci=T,type = "html")
```
# Generate plots








## Composite temperature dataset

```{r tab_mod2, include=T, results='asis'}
cat(htmlPreserve(tab2))
```

```{r temp_plot,include=TRUE}
#png(filename="composite_temp.png", width=180, height=100, units = "mm", res = 1200)
met_plot111 <- ggplot() +
  theme_bw() +
  geom_line(aes(x=ptime,y=avta,linetype=after),data = ta_composite) +
   labs(
    x = "Year", 
    y = "Daily average temperature (degC)",
    colour = "After installation"
    ) + theme(legend.position="none")
met_plot111
#dev.off()
```


## Whole house electrical demand profiles 2019

```{r dat2019, include=T}
v1_all <- dat2019_long %>% ggplot() +
  theme_bw()+
  geom_line(aes(x=Hours,y=value)) +
  facet_wrap(~Month.ct) +
  labs(
    x="Hour of day",
    y="Average electrical demand (kW)"
  )
v1_all
```


## Comparison of billing data with monitoring data


```{r billvsmonitor, include=TRUE}

#png(filename="step_consump.png", width=180, height=100, units = "mm", res = 1200)
plot2 <- ggplot() +
  theme_bw() +
  geom_step(data = monthly_sums,aes(x=ptime,y=cumEnergy_from_grid,linetype = Villa)) +
  geom_point(data = the_bill_data,aes(x=BillMonth,y=cumConsKWh,shape = Villa)) +
  scale_y_continuous(labels = label_comma()) +
   labs(
    x = "Year", 
    y = "Cumulative monthly consumption (kWh)"
    )
plot2
#dev.off()
```


## PV generation, demand and CDD

```{r plotting_ts,include=TRUE}
coeff <- 10

# keep every third string in a vector of strings 
everythird <- function(x){
  x <- sort(unique(x))
  x[seq(2, length(x), 3)] <- ""
  x[seq(3, length(x), 3)] <- ""
  x
}

plot3df <- joined_up_data_long %>% filter(ptime<"2024-06-01")

#png(filename="monthly_consump.png", width=180, height=100, units = "mm", res = 1200)
plot3 <- ggplot(plot3df) + 
  geom_col(aes(x=ym,y=value,fill = name),position = "stack") +
  theme(axis.text.x = element_text(angle = 90)) +
  facet_wrap(~Villa,labeller = top_string) +
  geom_path(aes(x=ym, y=cdd.adj*coeff),data = Composite_met_monthly_stats, group=1)+
  
  scale_y_continuous(
    
    # Features of the first axis
    name = "Monthly energy (kWh)",
    
    # Add a second axis and specify its features
    sec.axis = sec_axis(~./coeff, name="Monthly CDD")
  ) + 
  scale_x_discrete(labels = everythird(plot3df$ym), guide = guide_axis(angle = 90)) +
xlab("Year-month (YYYYMM)") +
  theme_bw()+
  #theme(axis.text.x = element_text(angle = 90, vjust = 1, hjust=1)) +
  scale_fill_discrete(name = "",labels=c("Total consumption", "PV generation"))

plot3
#dev.off()



```



## Variation of electrical demand with temperature

```{r elecDemand, include=TRUE}
#png(filename="cdd_fit.png", width=180, height=100, units = "mm", res = 1200)
met_plot22 <- composite_met_bill %>% ggplot() +
  theme_bw() +
  #geom_point(aes(x=cdd.adj,y=Villa_1,colour = after),data = composite_met_bill,shape = 1) + 
  geom_point(aes(x=cdd.adj,y=cons2,shape = after)) + 
  geom_line(aes(x=cdd.adj,y=mod_cons_l)) +
    labs(
    #x = bquote(paste("Cooling degree-days/month base ", 18.0 ~degree C,)), 
    x = bquote(paste("Cooling degree-days/month (base temperature ",18.0~degree,"C",")")),
    y = bquote(paste("Consumption ", (kWh/m^2/month)))
    ) +
#    scale_x_continuous(limits = c(0, 675)) +
#  scale_y_continuous(limits = c(0, 6000))+
  scale_shape(solid = FALSE,name = "",labels=c("Before PV", "After PV"))+
  facet_wrap(~Villa,labeller = top_string)
met_plot22
#dev.off()
```

```{r tab_mod1, include=T, results='asis'}
cat(htmlPreserve(tab1))
```

## PV performance

```{r PVperf,include=TRUE}
p2 <- the_data[the_data$ptime>="2023-06-03" & the_data$Villa==2 ,]$ifrac
p1 <- the_data[the_data$ptime< "2023-06-02" & the_data$Villa==2 ,]$ifrac
wilcox.test(p1,p2,alternative = "two.sided",paired=T)
ifracplot <- data.frame(p1,p2) %>% ggplot() + theme_bw() +
  geom_point(aes(x=p1,y=p2)) + 
  geom_abline(slope = 1,intercept = 0) +
    labs(
    x = "Generation/incident radiation Jun '22 - Jun '23 (kWh/kWh)",
    y = "Generation/incident radiation Jun '23 - Jun '24 (kWh/kWh)"
    ) +
    scale_x_continuous(limits = c(0, 0.2)) +
    scale_y_continuous(limits = c(0, 0.2))
ifracplot
```

# Financial analysis

## Parameters

```{r financeParams,include=TRUE,echo=TRUE}
discount_rate <- 0.15 # decimal
usd <- 3.75 # usd sar exchange rate 22 Nov 2024 (pegged to dollar)
tariff_exp <- 0.07/usd # current export tariff where applicable
tariff_low <- 0.18/usd # low rate import tariff (<6000kWh/month)
tariff_hi  <- 0.30/usd # high rate import tariff (>6000kWh/month)
tariff_eco <- 1.00/usd # theoretical 'economic' export tariff
tariff_vhi <- 1.50/usd # theoretical economic export tariff without other changes
oandm <- 150 # annual maintenance cost usd
cost_kWp <- 1250 # installed
capacity <- 15 # kWp
initial_capital <- cost_kWp * capacity
econ_scale <- 0.8 # 20% reduction in cost
tariff_increase <- 1.0 # i.e. NO increase
n_years <- 7 # 
```


## Consumption, generation and import

```{r consump,include=TRUE}
years_elapsed <- seq(1,n_years)

# average of 2 villas, 2 years
ann_consumption <- sum(the_data$Consumption)/4/1000 # kwh
ann_import      <- sum(the_data$Energy_from_grid)/4/1000 # kwh
ann_generation  <- sum(the_data$Consumed_directly)/4/1000 # kwh

# monthly means
ann_villa_month <- monthly_sums %>% filter(ptime < "2024-06-01") %>% group_by(mon=format(ptime, format="%m")) %>%
  summarize(Consumption=mean(Consumption),
            Generation=mean(Consumed_directly),
            Import=mean(Energy_from_grid), .groups = "drop"
           )

# use villa 2 for insolation as there are outliers on villa 1 irradiance data
mon_potential <- monthly_sums %>% filter(Villa==2) %>% filter(ptime < "2024-06-01") %>% group_by(mon=format(ptime, format="%m")) %>% summarize(Insolation_energy=mean(Insolation_energy),.groups="drop") # mean kWh array
ann_potential <- sum(mon_potential$Insolation_energy)
max_ifrac_row <- the_data %>% filter(Villa==2) %>% filter(ifrac==max(ifrac))
max_ifrac_row$ptime
max_ifrac_row$ifrac
ann_potential_pv <- ann_potential*max_ifrac_row$ifrac
ann_curtail <- ann_potential_pv - ann_generation
mon_potential %<>% mutate(potential_gen=Insolation_energy*max_ifrac_row$ifrac)

# overall villa consumption (using 2019 data aggregated to monthly)
whole_house <- data.frame(mon=sprintf("%02d",seq(1,12)),demand=c(3500,3400,4300,6000,6600,7000,8200,7700,7000,5900,3600,3400))
mon_dat <- bind_cols(mon_potential,ann_villa_month,whole_house)
mon_dat %<>% mutate(hi5M=if_else(demand>3000,demand-3000,0),
                    hi6M=if_else(demand>6000,demand-6000,0)
                    )
mon_dat %<>% mutate(hi5M_sav=if_else(hi5M>potential_gen,potential_gen,hi5M),
                    hi6M_sav=if_else(hi6M>potential_gen,potential_gen,hi6M)
                    )
ann_sav_6M <- sum(mon_dat$hi6M_sav)*(tariff_hi-tariff_low) # difference over threshold
ann_sav_5M <- sum(mon_dat$hi5M_sav)*(tariff_hi-tariff_low)
whole_house_import <- sum(mon_dat$demand) - sum(mon_dat$Generation)
```

## Cashflow, NPV and IRR 

```{r cashflow,echo=TRUE}
cashflow_bau     <- rep(                         ann_generation*tariff_low - oandm, n_years) # business as usual
cashflow_case1   <- rep(ann_curtail*tariff_exp + ann_generation*tariff_low - oandm, n_years) # Case I: export at current FiT
cashflow_case2   <- rep(ann_curtail*tariff_low + ann_generation*tariff_low - oandm, n_years) # Case II: export with FiT set at current import
cashflow_case3   <- rep(ann_curtail*tariff_hi  + ann_generation*tariff_low - oandm, n_years) # Case III: export with FiT set at current high rate
cashflow_case4   <- rep(ann_curtail*tariff_eco + tariff_increase*ann_generation*tariff_low - oandm, n_years) # Case IV: hypothetical FiT
cashflow_case5   <- rep(ann_curtail*tariff_eco + tariff_increase*(ann_generation*tariff_low + ann_sav_6M) - oandm, n_years) # Case V: hypothetical threshold applies at household level
cashflow_case6   <- rep(ann_curtail*tariff_eco + tariff_increase*(ann_generation*tariff_low + ann_sav_5M) - oandm, n_years) # Case VI: hypothetical lower threshold applies at household level
cashflow_hi      <- rep(ann_curtail*tariff_vhi + ann_generation*tariff_low - oandm, n_years) # FiT only
cashflow_zb      <- rep(                         ann_import*tariff_low     - oandm, n_years) # capped at zero bill for meter
cashflow_zb2     <- rep(                         whole_house_import*tariff_low - oandm, n_years) # capped at zero bill for household

npv_bau    <- NPV(cf0=initial_capital,            cf=cashflow_bau,    times=years_elapsed, i=discount_rate)
npv_case1  <- NPV(cf0=initial_capital,            cf=cashflow_case1,  times=years_elapsed, i=discount_rate)
npv_case2  <- NPV(cf0=econ_scale*initial_capital, cf=cashflow_case2,  times=years_elapsed, i=discount_rate)
npv_case3  <- NPV(cf0=econ_scale*initial_capital, cf=cashflow_case3,  times=years_elapsed, i=discount_rate)
npv_case4  <- NPV(cf0=econ_scale*initial_capital, cf=cashflow_case4,  times=years_elapsed, i=discount_rate)
npv_case5  <- NPV(cf0=econ_scale*initial_capital, cf=cashflow_case5,  times=years_elapsed, i=discount_rate)
npv_case6  <- NPV(cf0=econ_scale*initial_capital, cf=cashflow_case6,  times=years_elapsed, i=discount_rate)
npv_hi     <- NPV(cf0=initial_capital,            cf=cashflow_hi,     times=years_elapsed, i=discount_rate)
npv_zb     <- NPV(cf0=initial_capital,            cf=cashflow_zb,     times=years_elapsed, i=discount_rate)
npv_zb2    <- NPV(cf0=initial_capital,            cf=cashflow_zb2,    times=years_elapsed, i=discount_rate)

irr_bau    <- IRR(cf0=initial_capital,            cf=cashflow_case1,  times=years_elapsed)
irr_case1  <- IRR(cf0=initial_capital,            cf=cashflow_case2,  times=years_elapsed)
irr_case2  <- IRR(cf0=econ_scale*initial_capital, cf=cashflow_case3,  times=years_elapsed)
irr_case3  <- IRR(cf0=econ_scale*initial_capital, cf=cashflow_case3,  times=years_elapsed)
irr_case4  <- IRR(cf0=econ_scale*initial_capital, cf=cashflow_case4,  times=years_elapsed)
irr_case5  <- IRR(cf0=econ_scale*initial_capital, cf=cashflow_case5,  times=years_elapsed)
irr_case6  <- IRR(cf0=econ_scale*initial_capital, cf=cashflow_case6,  times=years_elapsed)
irr_hi     <- IRR(cf0=initial_capital,            cf=cashflow_hi,     times=years_elapsed)
irr_zb     <- IRR(cf0=initial_capital,            cf=cashflow_zb,     times=years_elapsed)
irr_zb2    <- IRR(cf0=initial_capital,            cf=cashflow_zb2,    times=years_elapsed)
```


```{r niceTable,include=TRUE}
print_irr <- function(x){
  if_else(length(x)==3,
          100*x[2],
          NULL)
}

fa <- data.frame(Case=c("BAU","I","II","III","IV","V","VI"),
                 NPV=c(npv_bau,npv_case1,npv_case2,npv_case3,npv_case4,npv_case5,npv_case6),
                 IRR=c(print_irr(irr_bau),print_irr(irr_case1),print_irr(irr_case2),
                       print_irr(irr_case3),print_irr(irr_case4),
                       print_irr(irr_case5),print_irr(irr_case6)))

kable(fa,digits = 2, format = "markdown") 


```





