johannes — Nov 23, 2012, 10:11 AM
#aids simulation
options(digits=2,scipen=10)
library(plyr)
library(MASS)
library(abind)
generateoutcomes=function(parameters,relshocksizeofsupply=1,correlationofshocks=0,nofshocks=100
){
shocked=mvrnorm(nofshocks,c(parameters[["alphad"]],parameters[["alphas"]]),
matrix(c(1,
correlationofshocks*relshocksizeofsupply^(1/2),
correlationofshocks*relshocksizeofsupply^(1/2),relshocksizeofsupply),nrow=2)
)
shocked=cbind(shocked,parameters[["demand"]])
shocked=cbind(shocked,parameters[["supply"]])
colnames(shocked)=c("alphad","alphas","demandPcoef","supplyQcoef")
solveshockedsystem=function(x){
#Q=AlphaD - demandPcoef*P --> AlphaD = Q+ demandPcoef*P
#P=AlphaS + supplyQcoef*Q --> AlphaS = -supplycoef*Q + P
a=matrix(c(1,x["demandPcoef"],-x["supplyQcoef"],1),nrow=2,byrow=TRUE)
solve(a,c(x["alphad"],x["alphas"])) #the first result is Q the second is P
}
equilibria=aaply(shocked,1,solveshockedsystem)
colnames(equilibria)=c("Q","P")
bb=cbind(equilibria,shocked[,1:2])
bb
}
#marketoutcomes=aaply(markets,1,generateoutcomes,relshocksizeofsupply=2,correlationofshocks=0,nofshocks=100)
evaluatebias=function(mo,parameters){
l2=lm(Q~P,data.frame(mo)) # estimate to compare -coef on P
estimate=l2$coefficients["P"]
bias=estimate - -parameters["demand"]
biasperc=(l2$coefficients["P"]/-parameters["demand"] -1)*100
res=c(bias,biasperc,estimate)
names(res)=c("bias absolute","bias percentage","estimate")
res
}
runmarket=function(market,relshocksizeofsupply,correlationofshocks,nofshocks){
mo=generateoutcomes(parameters=market,
relshocksizeofsupply=relshocksizeofsupply,
correlationofshocks=correlationofshocks,
nofshocks)
bias=evaluatebias(mo,parameters=market)
#plotmarket(mo,market)
bias
}
plotmarket=function(marketoutcome,market){
marketdata=marketoutcome[,1:2] #marketoutcomes
marketshocks=marketoutcome[,3:4] #market shocks alphad and alphas
rq=range(marketdata[,"Q"])
rp=range(marketdata[,"P"])
par(mar=c(5.1, 4.1, 4.1, 8.1), xpd=TRUE)
plot(rq,rp,type="n")
title(paste("market ",market["marketid"]))
points(marketdata[,"Q"],marketdata[,"P"])
#P= (alpha-Q)/suppcoef --> alpha/suppcoef - Q / suppcoef
#plot inverse demand and inverse supply curve
abline(market["alphad"]/market["demand"],-1/market["demand"],lty=3)
abline(market["alphas"],market["supply"],col="green",lty=3)
l=lm(Q~P,data.frame(marketdata)) #coef on price
estimate_alpha=l$coefficients[1]
estimate_coef=l$coefficients[2]
abline(estimate_alpha/-estimate_coef,1/estimate_coef,col="red")
bias=l$coefficients["P"] - -market["demand"]
bias
legend("topright",legend=c(paste("estimate d.cof",round(l$coefficients["P"],2)),
paste("dem coef bias:",round(bias,2)),
paste("demand coef",round(-market["demand"],2)),
paste("supply coef",round(market["supply"],2))),bg="transparent",ncol=1)
}
#############################################################################################################
inputs=list()
inputs$nofmarkets=200
inputs$demandcoefmean=.5 # coefficient on price of the demand side
inputs$demandcoefse=0
inputs$supplycoefmean=0 #coefficient on quantity of the inverse supply cruve
inputs$supplycoefse=0
attach(inputs)
#generate valid coefs between 0...~
demandcoefs= rnorm(nofmarkets+100,mean=mean(demandcoefmean),sd=demandcoefse)
demandcoefs=demandcoefs[demandcoefs >0 ][1:nofmarkets]
supplycoefs= rnorm(nofmarkets+100,mean=mean(supplycoefmean),sd=supplycoefse)
supplycoefs=supplycoefs[1:nofmarkets]
#for each pair of supply and demand solve system with N independent observations with shock
markets=as.matrix(data.frame(demand=demandcoefs,supply=supplycoefs,alphad=100,alphas=0,marketid=1:nofmarkets))
for (r in 1:10){
print("small sample 36 observations, show how bias in estimate varies with the relative size of supply shifts")
print(paste(r," times as big supply shift as demand shift"))
checkbias=aaply(markets,1,runmarket,relshocksizeofsupply=r,correlationofshocks=0,nofshocks=36)
colMeans(checkbias)
hist(checkbias[,"bias percentage"],main=paste("bias in% von Coef on P \n size of sup. shift",r))
}
[1] "small sample 36 observations, show how bias in estimate varies with the relative size of supply shifts"
[1] "1 times as big supply shift as demand shift"
[1] "small sample 36 observations, show how bias in estimate varies with the relative size of supply shifts"
[1] "2 times as big supply shift as demand shift"
[1] "small sample 36 observations, show how bias in estimate varies with the relative size of supply shifts"
[1] "3 times as big supply shift as demand shift"
[1] "small sample 36 observations, show how bias in estimate varies with the relative size of supply shifts"
[1] "4 times as big supply shift as demand shift"
[1] "small sample 36 observations, show how bias in estimate varies with the relative size of supply shifts"
[1] "5 times as big supply shift as demand shift"
[1] "small sample 36 observations, show how bias in estimate varies with the relative size of supply shifts"
[1] "6 times as big supply shift as demand shift"
[1] "small sample 36 observations, show how bias in estimate varies with the relative size of supply shifts"
[1] "7 times as big supply shift as demand shift"
[1] "small sample 36 observations, show how bias in estimate varies with the relative size of supply shifts"
[1] "8 times as big supply shift as demand shift"
[1] "small sample 36 observations, show how bias in estimate varies with the relative size of supply shifts"
[1] "9 times as big supply shift as demand shift"
[1] "small sample 36 observations, show how bias in estimate varies with the relative size of supply shifts"
[1] "10 times as big supply shift as demand shift"
# ####no case when supply shocks are as big as demand shocks
#
# checkbias=aaply(markets,1,runmarket,relshocksizeofsupply=1,correlationofshocks=0,nofshocks=100)
# colMeans(checkbias)
# hist(checkbias[,"bias percentage"],main="bias in% von Coef on P \n with non correlated shocks and \n equal as large supply shocks")
#
#
# checkbias=aaply(markets,1,runmarket,relshocksizeofsupply=10,correlationofshocks=0,nofshocks=100)
# colMeans(checkbias)
# hist(checkbias[,"bias percentage"])
#
# checkbias=aaply(markets,1,runmarket,relshocksizeofsupply=.1,correlationofshocks=0,nofshocks=100)
# colMeans(checkbias)
# hist(checkbias[,"bias percentage"])
#
#
# ##aslon as the supply curve is flat - the bias is very small, but with a steeper supply curve, the bias
# #is substantial
# inputs=list()
# inputs$nofmarkets=200
# inputs$demandcoefmean=.5 # coefficient on price of the demand side
# inputs$demandcoefse=0
# inputs$supplycoefmean=1 #coefficient on quantity of the inverse supply cruve
# inputs$supplycoefse=0
#
#
# attach(inputs)
# #generate valid coefs between 0...~
# demandcoefs= rnorm(nofmarkets+100,mean=mean(demandcoefmean),sd=demandcoefse)
# demandcoefs=demandcoefs[demandcoefs >0 ][1:nofmarkets]
# supplycoefs= rnorm(nofmarkets+100,mean=mean(supplycoefmean),sd=supplycoefse)
# supplycoefs=supplycoefs[1:nofmarkets]
# #for each pair of supply and demand solve system with N independent observations with shock
# markets=as.matrix(data.frame(demand=demandcoefs,supply=supplycoefs,alphad=0,alphas=0,marketid=1:nofmarkets))
#
#
# checkbias=aaply(markets,1,runmarket,relshocksizeofsupply=1,correlationofshocks=0,nofshocks=100)
# colMeans(checkbias)
# hist(checkbias[,"bias percentage"])
#
# # in this case of much bigger supply shifts, the bias is reduced
#
# checkbias=aaply(markets,1,runmarket,relshocksizeofsupply=16,correlationofshocks=0,nofshocks=100)
# colMeans(checkbias)
# hist(checkbias[,"bias percentage"])
#
# #if shocks are positively correlated in this case, the estimate looks like a good fit - but its way off the real value
# checkbias=aaply(markets,1,runmarket,relshocksizeofsupply=16,correlationofshocks=1,nofshocks=100)
# colMeans(checkbias)
# hist(checkbias[,"bias percentage"])
#
#
# ##what happens in the case with flat supply but correlated demand and supply?
# inputs=list()
# inputs$nofmarkets=100
# inputs$demandcoefmean=.5 # coefficient on price of the demand side
# inputs$demandcoefse=0
# inputs$supplycoefmean=0 #coefficient on quantity of the inverse supply cruve
# inputs$supplycoefse=0
#
#
# attach(inputs)
# #generate valid coefs between 0...~
# demandcoefs= rnorm(nofmarkets+100,mean=mean(demandcoefmean),sd=demandcoefse)
# demandcoefs=demandcoefs[demandcoefs >0 ][1:nofmarkets]
# supplycoefs= rnorm(nofmarkets+100,mean=mean(supplycoefmean),sd=supplycoefse)
# supplycoefs=supplycoefs[1:nofmarkets]
# #for each pair of supply and demand solve system with N independent observations with shock
# markets=as.matrix(data.frame(demand=demandcoefs,supply=supplycoefs,alphad=0,alphas=0,marketid=1:nofmarkets))
#
# #strong correlation
# checkbias=aaply(markets,1,runmarket,relshocksizeofsupply=1,correlationofshocks=.5,nofshocks=100)
# colMeans(checkbias)
# hist(checkbias[,"bias percentage"])
#
#
# #week correlation
# checkbias=aaply(markets,1,runmarket,relshocksizeofsupply=1,correlationofshocks=.1,nofshocks=100)
# colMeans(checkbias)
# hist(checkbias[,"bias percentage"])
#
# #negative correlation
# checkbias=aaply(markets,1,runmarket,relshocksizeofsupply=1,correlationofshocks=-.1,nofshocks=100)
# colMeans(checkbias)
# hist(checkbias[,"bias percentage"])
#
# #negative correlation -leads to a seizable bias, if shocks are of equal size
# checkbias=aaply(markets,1,runmarket,relshocksizeofsupply=2,correlationofshocks=-.5,nofshocks=100)
# colMeans(checkbias)
# hist(checkbias[,"bias percentage"])
#
#
# #if supply shocks are bigger and there is no correlation of the shocks - then the estimate is good
# checkbias=aaply(markets,1,runmarket,relshocksizeofsupply=2,correlationofshocks=0,nofshocks=100)
# colMeans(checkbias)
# hist(checkbias[,"bias percentage"])
#
#