library(lubridate)
##
## Attaching package: 'lubridate'
## The following objects are masked from 'package:base':
##
## date, intersect, setdiff, union
library(ggplot2)
library(lfe)
## Warning: package 'lfe' was built under R version 4.3.2
## Loading required package: Matrix
library(sf)
## Linking to GEOS 3.11.2, GDAL 3.6.2, PROJ 9.2.0; sf_use_s2() is TRUE
library(tmap)
## The legacy packages maptools, rgdal, and rgeos, underpinning the sp package,
## which was just loaded, will retire in October 2023.
## Please refer to R-spatial evolution reports for details, especially
## https://r-spatial.org/r/2023/05/15/evolution4.html.
## It may be desirable to make the sf package available;
## package maintainers should consider adding sf to Suggests:.
## The sp package is now running under evolution status 2
## (status 2 uses the sf package in place of rgdal)
library(data.table)
##
## Attaching package: 'data.table'
## The following objects are masked from 'package:lubridate':
##
## hour, isoweek, mday, minute, month, quarter, second, wday, week,
## yday, year
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:data.table':
##
## between, first, last
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats 1.0.0 ✔ stringr 1.5.0
## ✔ purrr 1.0.1 ✔ tibble 3.2.1
## ✔ readr 2.1.4 ✔ tidyr 1.3.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::between() masks data.table::between()
## ✖ tidyr::expand() masks Matrix::expand()
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::first() masks data.table::first()
## ✖ data.table::hour() masks lubridate::hour()
## ✖ data.table::isoweek() masks lubridate::isoweek()
## ✖ dplyr::lag() masks stats::lag()
## ✖ dplyr::last() masks data.table::last()
## ✖ data.table::mday() masks lubridate::mday()
## ✖ data.table::minute() masks lubridate::minute()
## ✖ data.table::month() masks lubridate::month()
## ✖ tidyr::pack() masks Matrix::pack()
## ✖ data.table::quarter() masks lubridate::quarter()
## ✖ data.table::second() masks lubridate::second()
## ✖ purrr::transpose() masks data.table::transpose()
## ✖ tidyr::unpack() masks Matrix::unpack()
## ✖ data.table::wday() masks lubridate::wday()
## ✖ data.table::week() masks lubridate::week()
## ✖ data.table::yday() masks lubridate::yday()
## ✖ data.table::year() masks lubridate::year()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(zoo)
##
## Attaching package: 'zoo'
##
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
library(slider)
## Warning: package 'slider' was built under R version 4.3.2
library(timeDate)
##### Robust IV Estimation for Spend with Air temperature directly (rather than Extreme Heat)#####
################################################################################################
################################################################################################
filesave <-paste0("E:/POI_gridmet_daily_weather_data/7cities_heat_spend_visit_climate",".RData")
load(filesave)
Q1 <- quantile(df$spend_pertransaction, 0.25, na.rm = TRUE) # 1st Quartile
Q3 <- quantile(df$spend_pertransaction, 0.75, na.rm = TRUE) # 3rd Quartile
IQR <- Q3 - Q1 # Interquartile Range
# Define outlier thresholds
lower_bound <- Q1 - 1.5 * IQR
upper_bound <- Q3 + 1.5 * IQR
df <- df[df$spend_pertransaction >= lower_bound & df$spend_pertransaction <= upper_bound, ]
#IV estimation
files <- list.files(path = "E:/POI_gridmet_daily_weather_data/grass", full.names = TRUE) #read grass evaporation data
df_grass_eva<-data.frame()
for (i in 1:40){
load(files[i])
print(i)
daily_grass$eva <- daily_grass$grass$`pet_2019-05-01`
daily_grass$grass <-NULL
current_df <- merge(df, daily_grass, by=c("PLACEKEY","date"))
df_grass_eva <-bind_rows(df_grass_eva,current_df)
}
## [1] 1
## [1] 2
## [1] 3
## [1] 4
## [1] 5
## [1] 6
## [1] 7
## [1] 8
## [1] 9
## [1] 10
## [1] 11
## [1] 12
## [1] 13
## [1] 14
## [1] 15
## [1] 16
## [1] 17
## [1] 18
## [1] 19
## [1] 20
## [1] 21
## [1] 22
## [1] 23
## [1] 24
## [1] 25
## [1] 26
## [1] 27
## [1] 28
## [1] 29
## [1] 30
## [1] 31
## [1] 32
## [1] 33
## [1] 34
## [1] 35
## [1] 36
## [1] 37
## [1] 38
## [1] 39
## [1] 40
########stage I########
first_stage <- felm(tmax ~eva+ pr+rmax + wind|
PLACEKEY + date|0,
data = df_grass_eva)
summary(first_stage)
##
## Call:
## felm(formula = tmax ~ eva + pr + rmax + wind | PLACEKEY + date | 0, data = df_grass_eva)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.0509 -0.8611 -0.0002 0.9412 9.4252
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## eva 1.6739073 0.0017922 934.0 <2e-16 ***
## pr 0.0214800 0.0001628 131.9 <2e-16 ***
## rmax -0.0085049 0.0001636 -52.0 <2e-16 ***
## wind -0.6617640 0.0016792 -394.1 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.652 on 2163342 degrees of freedom
## Multiple R-squared(full model): 0.8728 Adjusted R-squared: 0.8722
## Multiple R-squared(proj model): 0.4571 Adjusted R-squared: 0.4547
## F-statistic(full model): 1572 on 9441 and 2163342 DF, p-value: < 2.2e-16
## F-statistic(proj model): 4.554e+05 on 4 and 2163342 DF, p-value: < 2.2e-16
##get the fitted values
fixed_effects <- getfe(first_stage)$fe
model_coefs <- coef(first_stage) # Extracting model coefficients
model_mat <- model.matrix(~ 0+eva + pr + rmax + wind, data = df_grass_eva) #create matrix
linear_predictor <- model_mat %*% model_coefs
df_grass_eva$fitted_heat <-linear_predictor
#########stage II########
second_stage <- felm(spend_pertransaction ~fitted_heat+ pr+rmax + wind|
PLACEKEY + date|0,
data = df_grass_eva)
summary(second_stage)
##
## Call:
## felm(formula = spend_pertransaction ~ fitted_heat + pr + rmax + wind | PLACEKEY + date | 0, data = df_grass_eva)
##
## Residuals:
## Min 1Q Median 3Q Max
## -37.231 -5.267 -1.515 3.529 48.543
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## fitted_heat -0.2489432 0.0059550 -41.804 <2e-16 ***
## pr -0.0012748 0.0008868 -1.438 0.151
## rmax 0.0236380 0.0009446 25.025 <2e-16 ***
## wind 0.0670961 0.0077421 8.666 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 9.19 on 2163342 degrees of freedom
## Multiple R-squared(full model): 0.3186 Adjusted R-squared: 0.3156
## Multiple R-squared(proj model): 0.003584 Adjusted R-squared: -0.0007649
## F-statistic(full model):107.1 on 9441 and 2163342 DF, p-value: < 2.2e-16
## F-statistic(proj model): 1945 on 4 and 2163342 DF, p-value: < 2.2e-16
##### Robust IV Estimation for Visit with Air temperature directly (rather than Extreme Heat)#####
##################################################################################################
##################################################################################################
#### filter out outliers
filesave <-paste0("E:/POI_gridmet_daily_weather_data/7cities_heat_spend_visit_climate",".RData")
load(filesave)
Q1 <- quantile(df$visit, 0.25, na.rm = TRUE) # 1st Quartile
Q3 <- quantile(df$visit, 0.75, na.rm = TRUE) # 3rd Quartile
IQR <- Q3 - Q1 # Interquartile Range
lower_bound <- Q1 - 1.5 * IQR
upper_bound <- Q3 + 1.5 * IQR
df <- df[df$visit >= lower_bound & df$visit <= upper_bound, ]
df <- df[!is.na(df$NAME), ]
#IV estimation
files <- list.files(path = "E:/POI_gridmet_daily_weather_data/grass", full.names = TRUE) #read grass evaporation data
df_grass_eva<-data.frame()
for (i in 1:40){
load(files[i])
print(i)
daily_grass$eva <- daily_grass$grass$`pet_2019-05-01`
daily_grass$grass <-NULL
current_df <- merge(df, daily_grass, by=c("PLACEKEY","date"))
df_grass_eva <-bind_rows(df_grass_eva,current_df)
}
## [1] 1
## [1] 2
## [1] 3
## [1] 4
## [1] 5
## [1] 6
## [1] 7
## [1] 8
## [1] 9
## [1] 10
## [1] 11
## [1] 12
## [1] 13
## [1] 14
## [1] 15
## [1] 16
## [1] 17
## [1] 18
## [1] 19
## [1] 20
## [1] 21
## [1] 22
## [1] 23
## [1] 24
## [1] 25
## [1] 26
## [1] 27
## [1] 28
## [1] 29
## [1] 30
## [1] 31
## [1] 32
## [1] 33
## [1] 34
## [1] 35
## [1] 36
## [1] 37
## [1] 38
## [1] 39
## [1] 40
########stage I########
first_stage <- felm(tmax ~eva+ pr+rmax + wind|
PLACEKEY + date|0,
data = df_grass_eva)
summary(first_stage)
##
## Call:
## felm(formula = tmax ~ eva + pr + rmax + wind | PLACEKEY + date | 0, data = df_grass_eva)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.8316 -0.8628 0.0014 0.9446 9.4553
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## eva 1.6779341 0.0018309 916.44 <2e-16 ***
## pr 0.0212695 0.0001670 127.35 <2e-16 ***
## rmax -0.0082312 0.0001662 -49.52 <2e-16 ***
## wind -0.6593237 0.0017180 -383.77 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.649 on 2069307 degrees of freedom
## Multiple R-squared(full model): 0.8729 Adjusted R-squared: 0.8723
## Multiple R-squared(proj model): 0.4592 Adjusted R-squared: 0.4568
## F-statistic(full model): 1595 on 8905 and 2069307 DF, p-value: < 2.2e-16
## F-statistic(proj model): 4.392e+05 on 4 and 2069307 DF, p-value: < 2.2e-16
##get the fitted values
fixed_effects <- getfe(first_stage)$fe
model_coefs <- coef(first_stage) # Extracting model coefficients
model_mat <- model.matrix(~ 0+eva + pr + rmax + wind, data = df_grass_eva) #create matrix
linear_predictor <- model_mat %*% model_coefs
df_grass_eva$fitted_heat <-linear_predictor
#########stage II########
second_stage <- felm(visit ~fitted_heat+ pr+rmax + wind|
PLACEKEY + date|0,
data = df_grass_eva)
summary(second_stage)
##
## Call:
## felm(formula = visit ~ fitted_heat + pr + rmax + wind | PLACEKEY + date | 0, data = df_grass_eva)
##
## Residuals:
## Min 1Q Median 3Q Max
## -37.097 -3.042 -0.382 2.429 44.157
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## fitted_heat -0.3413138 0.0039262 -86.931 < 2e-16 ***
## pr -0.0031763 0.0005886 -5.396 6.8e-08 ***
## rmax 0.0139111 0.0006204 22.423 < 2e-16 ***
## wind -0.1040032 0.0051229 -20.302 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.932 on 2069307 degrees of freedom
## Multiple R-squared(full model): 0.6491 Adjusted R-squared: 0.6476
## Multiple R-squared(proj model): 0.01131 Adjusted R-squared: 0.007052
## F-statistic(full model):429.9 on 8905 and 2069307 DF, p-value: < 2.2e-16
## F-statistic(proj model): 5916 on 4 and 2069307 DF, p-value: < 2.2e-16