#########################nulstiller miljø###################

rm(list=ls())

###"########################### kalder libraries   ##############################

pack<-c("car","sandwich","lmtest","RColorBrewer","mgcv","foreign","xtable"
        ,"AER","stargazer", "MASS","jpeg")

lapply(pack, require, character.only=T)

rm(pack)



############################################# setting working path##########################
data_path<-"C:/Users/Schweppes/Desktop/Econometrics"
result_path<-"C:/Users/Schweppes/Desktop/Econometrics"
graphs_path<-"C:/Users/Schweppes/Desktop/Econometrics"

################# reading in adata########################################################


setwd(data_path)
dir()

load("adata")
desc

################################################################################################


###########################################################
### Tilføjer grid og tilhørende variable
###########################################################


adata$tilex = ((adata$x - min(adata$x)) %/% (((max(adata$x) + 1) - min(adata$x)) / 10)) * 1000 + 1000
adata$tiley = ((adata$y - min(adata$y)) %/% (((max(adata$y) + 1) - min(adata$y)) / 10)) * 10 + 10
adata$grid = as.factor(adata$tilex + adata$tiley)

adata$in12 = adata$in1 + adata$in2

max(adata$x)
min(adata$y)
options(scipen=3)
###########################################################

mean(adata$price)

UndersogLog<-lm(log(price)~area
                +age
                +bathrooms
                +rebuild90
                +rebuild00
                +roof_cemen
                +roof_fiber
                +roof_board
                +roof_flat
                +roof_tile
                +railways
                +I(railways^2)
                +big_roads
                +I(big_roads^2)
                +in12
                +in3
                +in4
                +allgoods_WA_1000m
                +divers_shopservice_WA_1000m
                +grid
                ,adata)

Undersog<-lm(price~area
             +age
             +bathrooms
             +rebuild90
             +rebuild00
             +roof_cemen
             +roof_fiber
             +roof_board
             +roof_flat
             +roof_tile
             +railways
             +I(railways^2)
             +big_roads
             +I(big_roads^2)
             +in12
             +in3
             +in4
             +allgoods_WA_1000m
             +divers_shopservice_WA_1000m
             +grid
             ,adata)
summary(Undersog)


bptest(Undersog)
coeftest(Undersog, vcov = vcovHC(Undersog, type = "HC"))
#waldtest(Undersog,UndersogNy, vcov = vcovHC)
resettest(Undersog)

par(mfrow=c(2,2))
hist(UndersogLog$residuals)
hist(Undersog$residuals)
qqnorm(residuals(UndersogLog))
qqline(residuals(UndersogLog))
qqnorm(residuals(Undersog))
qqline(residuals(Undersog))


length(adata$price[adata$in1==1])

nytHus1 <- data.frame(area=110,rebuild80=0,rebuild90=0,rebuild00=1,roof_cemen=1
                     ,roof_fiber=0,roof_board=0,roof_flat=0,roof_tile=0,grid="6050"
                     ,in1=0,in12=0,in3=0,in4=1,allgoods_WA_1000m=1, railways = 1120, big_roads = 280
                     ,divers_shopservice_WA_1000m=1,age=1990,bathrooms=2,floor=mean(adata$floor),date_num=mean(adata$date_num))

nytHus2 <- data.frame(area=110,rebuild80=0,rebuild90=0,rebuild00=1,roof_cemen=1
                      ,roof_fiber=0,roof_board=0,roof_flat=0,roof_tile=0,grid="10050"
                      ,in1=0,in12=0,in3=0,in4=1,allgoods_WA_1000m=1, railways = 1120, big_roads = 280
                      ,divers_shopservice_WA_1000m=1,age=1990,bathrooms=2,floor=mean(adata$floor),date_num=mean(adata$date_num))

predict(Undersog,nytHus1,interval="predict")
predict(Undersog,nytHus2,interval="predict")

exp((mean(Undersoglog$res^2))/2)*exp(predict(Undersoglog,nytHus1))

####################Heatmap for samlet datasæt###########################
par(mfrow=c(2,2))
library(plot3D)
neo = matrix(nrow = 10,ncol = 10)
for (i in c(1:10)){
  for(l in c(1:10)){
    neo[i,l] = if (length(adata$price[adata$tilex == (l*1000) & adata$tiley == (110-10*(i))]) == 0) {NA} else {length(adata$price[adata$tilex == (l*1000) & adata$tiley == (110-10*(i))])}
  }
}
neo = t(apply(neo, 2, rev)) ###roterer heatmap
image2D(neo, border = "black")
title(main="densitet over samlet datasæt", col.main="red")


###############################################################################

####################Heatmap for % af ældre mennesker###########################
library(plot3D)
neo = matrix(nrow = 10,ncol = 10)
for (i in c(1:10)){
  for(l in c(1:10)){
    neo[i,l] = length(adata$price[adata$AgeMin61 == 1  & adata$tilex == (l*1000) & adata$tiley == (110-10*(i))]) / length(adata$price[adata$tilex == (l*1000) & adata$tiley == (110-10*(i))])
  }
}
neo = t(apply(neo, 2, rev)) ###roterer heatmap
image2D(neo, border = "black")
title(main="%-del af folk over 60", col.main="red")
###############################################################################

####################Heatmap for gennemsnitlig pris###########################
library(plot3D)
neo = matrix(nrow = 10,ncol = 10)
for (i in c(1:10)){
  for(l in c(1:10)){
    neo[i,l] = mean(adata$price[adata$tilex == (l*1000) & adata$tiley == (110-10*(i)) & length(adata$price[adata$tilex == (l*1000) & adata$tiley == (110-10*(i))]) >=5])
  }
}
neo = t(apply(neo, 2, rev)) ###roterer heatmap
image2D(neo, border = "black")
title(main="gnm-snit af pris", col.main="red")
###############################################################################
v = c()


for (i in c(1:10)){
  for(l in c(1:10)){
    v = c(v,neo[i,l])
  }
}
max(v,na.rm =TRUE)
max(v,na.rm =TRUE)-min(v,na.rm =TRUE)
min(v,na.rm =TRUE)

####################Heatmap for gennemsnitlig liggetid###########################
library(plot3D)
neo = matrix(nrow = 10,ncol = 10)
for (i in c(1:10)){
  for(l in c(1:10)){
    neo[i,l] = mean(adata$date_num[adata$tilex == (l*1000) & adata$tiley == (110-10*(i))])
  }
}
neo = t(apply(neo, 2, rev)) ###roterer heatmap
image2D(neo, border = "black")
title(main="gnm-snit af liggetid i dage", col.main="red")
###############################################################################


#set.seed(1)
#m <- cbind( A=sample(Undersog$residuals, 1000, replace=FALSE), B=sample(Undersog$residuals, 1000, replace=FALSE) )
#colMeans(m)

#t.test(m[,1],m[,2], var.equal=TRUE, paired=FALSE)


###########Tester korrelation##############################
testvif = vif(Undersog)

testvif^2

adatacor = subset(adata, select = c(price
                                    ,area
                                    ,age
                                    ,bathrooms
                                    ,rebuild90
                                    ,rebuild00
                                    ,roof_cemen
                                    ,roof_fiber
                                    ,roof_board
                                    ,roof_flat
                                    ,roof_tile
                                    ,railways
                                   #,I(railways^2)
                                    ,big_roads
                                    #,I(big_roads^2)
                                    ,in12
                                    ,in3
                                    ,in4
                                    ,allgoods_WA_1000m
                                    ,divers_shopservice_WA_1000m))
res = cor(adatacor)
round(res,2)
##############################################################