#HANDS-ON INTERMEDIATE ECONOMETRICS USING R 
# Second Edition, 2022, World Scientific Publishers
#by Hrishikesh D Vinod 
#Following snippets are from Chapter 4
#If you use them, please give credit to Prof. Vinod



#R4.2.1  Lorenz Curve Plots from data Using R
library(ineq)
# 1998 USA income distribution in 10 classes
# x=vector of class means, n=vector of class frequencies
x <- c(7, 18, 32, 58)
n <- c(20, 20, 20, 20)
# compute the Lorenz curve
Lc1 <- Lc(x, n=n)
plot(Lc1,main="Lorenz curve (solid line), lognormal (dashed)")#quebrada
# add the theoretic Lorenz curve of a Lognormal-
#distribution with variance 0.78
lines(Lc.lognorm, parameter=0.78)#suavizada
Gini(x) # 0.3630435


#R4.2.2  Gini, Atkinson, Entropy, etc. coeff Using R
library(ineq)
# data vector of incomes
x <- c(641, 1463, 2445, 3438, 5401, 6392, 8304, 9004)
# compute Gini coefficient
ineq(x); Gini(x)
# compute Atkinson coefficient with parameter=0.5
ineq(x, parameter=0.5, type="Atkinson")
RS(x)
Theil(x, parameter = 0)
Kolm(x, parameter = 1)
var.coeff(x, square = FALSE)
entropy(x, parameter = 0.5)



#R4.2.3 general R function for Prelec Transformation
prelec=function(n, alp)
  #inputs: size n & value of alpha between 0 and 1
{x=1:n; p=rep(0,length(x)); wp=p;
#loop to calculate cumulative probabilities
for (i in 1:n) {p[i]=x[i]/n; lnp=log(p[i])
wp[i]=exp(-(-lnp)^alp)}
wpdif=wp # to get the first value of wpdif[1] as wp[1]
for (i in 2:n) {
  wpdif[i]=wp[i]-wp[i-1]} #compute first differences
list(x=x,p=p,wp=wp,wpdif=wpdif)  }



#R4.3.1  Beta Density Illustrates first-order
#stochastic dominance Using R
x=seq(0,1,by=0.01)# define the range of x
plot(x, dbeta(x,2,5), typ="l")
# add the dashed line
lines(x, dbeta(x,3,6), typ="l", lty=2)


#R4.3.2 R function for matrices I sub f and I sub big F
ibothf=function(dj){
  n=length(dj); ibf=diag(dj)
  isf=diag(n); ibf[1,1]=dj[1]
  for (i in 2:n){for (j in 1: (i-1)) {
    ibf[i,j]=dj[j]+dj[j+1]; isf[i,j]=1 }}
  list(I.smallf=isf, I.bigf=0.5*ibf)}



#R4.3.3 sort x matrix x by column j
sort.matrix =function(x,j){
  y=x[sort.list(x[,j]),]
  return(y)  }



#R4.3.4 general R function for weighted pa & pb & dj
#must have #R4.3.3 & #R4.2.3 in memory for this to work
wtdpapb=function(xa,xb, alp=1){
  #input: alpha and excess returns for mutual fund A & B
  #output is weighted pa and pb vectors (wpa, wpb) and dj vector
  Ta=length(xa); Tb=length(xb);k=Ta+Tb
  pa0=rep(1/Ta, Ta); pb0=rep(1/Tb, Tb)
  xpapb=matrix(0,k,3)
  for (i in 1:Ta){xpapb[i,1]=xa[i]
  xpapb[i,2]=pa0[i]}
  for (i in 1:Tb){xpapb[Ta+i,1]=xb[i]
  xpapb[Ta+i,3]=pb0[i]}
  pra=prelec(n=Ta+Tb, alp)
  prb=prelec(n=Ta+Tb, alp)
  sm=sort.matrix(xpapb,1) #sort on first col of excess returns
  pa=sm[,2]; pb=sm[,3]
  #print("lengths of pa and pra$wpdif")
  #print(c(length(pa), length(pra$wpdif)))
  wpa=pa*pra$wpdif; wpb=pb*prb$wpdif
  #weighted pa use first differences of prelec weights
  dj=sm[,1]-sm[1,1]  #deviations from the minimum
  #at sm[1,1] (sm is sorted)
  list(wpa=wpa, wpb=wpb, dj=dj)}



#R4.3.5 general function for 4 stochastic dominance orders
stochdom =function( wpa, wpb, I.smallf, I.bigf) {
  # input weighted pa and pb  IsubF and I sub f matrices
  if (length(wpa) != nrow(I.smallf)) print("wrong rows")
  if (length(wpa) != nrow(I.bigf)) print("wrong rows I sub F")
  sd1=I.bigf %*% I.smallf %*% (wpa-wpb)
  sd2=I.bigf %*% sd1
  sd3=I.bigf %*% sd2
  sd4=I.bigf %*% sd3
  list( sd1=sd1,sd2=sd2, sd3=sd3,sd4=sd4) }



#R4.3.6 general function for portfolio comparisons
#Place #R4.3.2 to #R4.3.5 in R memory before this.
comp.portfo = function(xa, xb, alp){
  if (alp<0.0001){ alp=0.0001; print(" alp reset at 0.0001")}
  #Be sure to remove all missing data
  # or NAs from the portfolio excess return data
  xaa=na.omit(xa)
  xbb=na.omit(xb)
  gp=wtdpapb(xaa,xbb, alp)  #function above
  ibo=ibothf(gp$dj)
  stdo=stochdom(gp$wpa, gp$wpb, ibo$I.bigf, ibo$I.smallf)
  out=cbind(stdo$sd1,stdo$sd2, stdo$sd3, stdo$sd4)
  #print(out)
  #column sums for 4 orders of stochastic dominance
  ousum=apply(out,2,sum)
  return(ousum)}   #  #  Function ends
#  #  #  Simple example  begins
alp=0.5
set.seed(45); xa=runif(10,2,5)
xb=runif(10,2.2,5.1); matplot(cbind(xa,xb))
comp.portfo(xa,xb,alp=.5)
# 1.403913  10.435054  56.622903 252.967237