#
#  Average derivatives (1997-2006) of the factors in Barigozzi & Moneta "Identifying the independent sources of consumption variation"
#
#  Results are inserted in Tables 9 of the paper.
#
# You need to put the three identifyied factor in a csv format.
# for example the factors which refer to the analysis with 2-4 members of households are saved in a csv file called "factors_24member.csv"
# This file is simply a three-columns matrix where each column is a factor. The columns are named "f1", "f2", and "f3"
#

rm(list = ls())
library(sm)
library(KernSmooth)
#source("C:\\Documents and Settings\\User\\Documenti\\factors\\Rcodes\\funktionen.r")
#fs<-read.table("U:\\factors\\Rcodes\\factors.txt")   # old data used in the previous version submitted to RES
source("/home/alessio/Dropbox/factors/Rcodes/funktionen.r")

fs<-read.csv("/home/alessio/Dropbox/BM_Engel_curves/prog/r4_revision/factors_24members.csv",      # <- <- check the year used!!!
             dec=".",sep = ",", header = FALSE, col.names=c("f1", "f2", "f3")) 
te<-1:100
E1<-cbind(te, fs[,1])    # first factor
E2<-cbind(te, fs[,2])    # second factor
E3<-cbind(te, fs[,3])    # third factor
### anylsis average derivative with Haerdle-Stoker method
D1<-fun_avderiv(E1)
D2<-fun_avderiv(E2)
D3<-fun_avderiv(E3)


### standard deviations

SDV1<- fun_avd_se(E1, niter=1000)
SDV2<- fun_avd_se(E2, niter=1000)
SDV3<- fun_avd_se(E3, niter=1000)
# Wald test
W1<- D1^2 / SDV1^2
W2<- D2^2 / SDV2^2
W3<- D3^2 / SDV3^2

### the following analysis is made for two quantiles ("poor" vs. "rich")
E1_1<-E1[1:50,]
E1_2<-E1[51:100,]

E2_1<-E2[1:50,]
E2_2<-E2[51:100,]

E3_1<-E3[1:50,]
E3_2<-E3[51:100,]

D1_1<-fun_avderiv(E1_1)
D1_2<-fun_avderiv(E1_2)

D2_1<-fun_avderiv(E2_1)
D2_2<-fun_avderiv(E2_2)

D3_1<-fun_avderiv(E3_1)
D3_2<-fun_avderiv(E3_2)

SDV1_1<- fun_avd_se(E1_1, niter=1000)
SDV1_2<- fun_avd_se(E1_2, niter=1000)

SDV2_1<- fun_avd_se(E2_1, niter=1000)
SDV2_2<- fun_avd_se(E2_2, niter=1000)

SDV3_1<- fun_avd_se(E3_1, niter=1000)
SDV3_2<- fun_avd_se(E3_2, niter=1000)


W1_1<- D1_1^2 / SDV1_1^2
W1_2<- D1_2^2 / SDV1_2^2

W2_1<- D2_1^2 / SDV2_1^2
W2_2<- D2_2^2 / SDV2_2^2

W3_1<- D3_1^2 / SDV3_1^2
W3_2<- D3_2^2 / SDV3_2^2

#calcolo p-values
p1 <-  1 - pchisq(W1, df=1)
p2 <-  1 - pchisq(W2, df=1)
p3 <-  1 - pchisq(W3, df=1)

p1_1<-  1 - pchisq(W1_1, df=1)
p1_2<-  1 - pchisq(W1_2, df=1)

p2_1<-  1 - pchisq(W2_1, df=1)
p2_2<-  1 - pchisq(W2_2, df=1)

p3_1<-  1 - pchisq(W3_1, df=1)
p3_2<-  1 - pchisq(W3_2, df=1)



# list of results
D1
D2
D3

D1_1
D1_2

D2_1
D2_2

D3_1
D3_2

W1
W2
W3

W1_1
W1_2

W2_1
W2_2

W3_1
W3_2

#p-values
p1
p2
p3

p1_1
p1_2

p2_1
p2_2

p3_1
p3_2
