Modelado de Precios con la base de Properati

Laboratorio de Datos

Carga de Librerías, Lectura del dataset y su estructura

# leo el archivo ar_properties 
library(tidyverse) # entorno tidy
library(dplyr) # manejo de datos
library(GGally) # scatterplots multiples
library(rgl) # para graficos 3D
datos1a <- read_csv("/home/andresfaral/Dropbox/Labo de Datos/ar_properties.csv") # Acá completen con su propio PATH al archivo

── Column specification ──────────────────────────────────────────────────────────────────────
cols(
  .default = col_character(),
  start_date = col_date(format = ""),
  end_date = col_date(format = ""),
  created_on = col_date(format = ""),
  lat = col_double(),
  lon = col_double(),
  l6 = col_logical(),
  rooms = col_double(),
  bedrooms = col_double(),
  bathrooms = col_double(),
  surface_total = col_double(),
  surface_covered = col_double(),
  price = col_double()
)
ℹ Use `spec()` for the full column specifications.

datos1a

Aplicando filtros

datos1d <- datos1a %>% 
                   # Me quedo con los que pertenecen a Argentina, Capital Federal y Boedo
            filter(l1 == "Argentina", 
                   l2 == "Capital Federal",
                   l3=="Boedo",
                   # cuyo precio este en dolares 
                   currency == "USD", 
                   # propiedad tipo Casa
                   property_type %in% c("Casa"),
                   # operaciones de venta
                   operation_type == "Venta") %>% 
            dplyr::select(id, l3, surface_total, surface_covered, price) %>% mutate(Precio=price,Sup=surface_covered,Fondo=surface_total-surface_covered) %>% dplyr::select(Sup,Fondo,Precio) %>%  filter(Fondo>=0) %>% na.omit()
datos1d

NA

Modelado de Datos

# Boedo
summary(datos1d)

      Sup            Fondo            Precio      
 Min.   : 35.0   Min.   :  0.00   Min.   : 90000  
 1st Qu.:136.5   1st Qu.:  7.00   1st Qu.:261250  
 Median :190.0   Median : 37.00   Median :300000  
 Mean   :190.1   Mean   : 85.47   Mean   :347997  
 3rd Qu.:251.0   3rd Qu.:180.00   3rd Qu.:450000  
 Max.   :332.0   Max.   :368.00   Max.   :590000  

ggpairs(datos1d)

# ajuste
ajusM2<-lm(Precio~Sup+Fondo,data=datos1d) # modelo lineal multiple
coe<-coef(ajusM2) # coeficientes
coe

(Intercept)         Sup       Fondo 
133400.9620    863.4824    590.5066 

summary(ajusM2)

Call:
lm(formula = Precio ~ Sup + Fondo, data = datos1d)

Residuals:
    Min      1Q  Median      3Q     Max 
-174199  -44261   11610   41648  230873 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 133401.0    27584.3   4.836 8.47e-06 ***
Sup            863.5      152.3   5.671 3.53e-07 ***
Fondo          590.5      108.7   5.434 8.85e-07 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 78620 on 65 degrees of freedom
Multiple R-squared:  0.6426,    Adjusted R-squared:  0.6316 
F-statistic: 58.43 on 2 and 65 DF,  p-value: 3.01e-15

plot3d(ajusM2,size=15,col="blue")
# Defino las variables
Sup<-datos1d$Sup
Fondo<-datos1d$Fondo
Precio<-datos1d$Precio

grafico bueno

predichos2<-fitted.values(ajusM2)
coefs <- coef(ajusM2)
a <- coefs["Sup"]
b <- coefs["Fondo"]
cc <- -1
d<- coefs["(Intercept)"]

par3d(windowRect = c(0, 0, 800, 800)) # make the window large
par3d(zoom = 1.1) # larger values make the image smaller
#plot3d(datos,col=colores[clases],size=2)
plot3d(Sup,Fondo,Precio, 
       type="s", size=1,col="red",pch="16", xlab="Sup", 
       ylab="Fondo", zlab="Precio")
planes3d(a, b, cc, d,  col = 'red', alpha = 0.2)
rgl.snapshot('/home/andresfaral/Dropbox/Labo de Datos/foto1.png')
plot3d(Sup,Fondo,predichos2, type="s", size=1,
       col="pink",pch="16",  xlab="Sup", 
       ylab="Fondo", zlab="Precio",add=T)
rgl.snapshot('/home/andresfaral/Dropbox/Labo de Datos/foto2.png')
segments3d(x=as.vector(rbind(Sup,Sup)),y=as.vector(rbind(Fondo,Fondo)),z=as.vector(rbind(Precio,predichos2)),col="darkred")
rgl.snapshot('/home/andresfaral/Dropbox/Labo de Datos/foto3.png')

Funcion de Evaluacion

# Perdida cuadratica
Eval <- function(mu, alfa, beta) {
  salida<-mean((Precio-mu-alfa*Sup-beta*Fondo)^2)
  return(t(salida))
}
Eval(100000,2000,1000)

            [,1]
[1,] 64555965000

Optimizacion Aleatoria

facred<-0.9999 # factor de reduccion de la ventana
facred.acu<-1 # factor de reduccion acumulado
toler<-1/1e5 # umbral de tolerancia
# rangos
rango.mu<- 50000 # rango inicial de mu
rango.alfa<- 500 # rango inicial de alfa
rango.beta<- 500 # rango inicial de beta
# parametros iniciales
mu<- 100000 # valor inicial de mu
alfa<- 1000 # valor inicial de alfa
beta<- 1000 # valor inicial de beta
# parametros mejores
mejor.mu<-mu
mejor.alfa<-alfa
mejor.beta<-beta
mejor.eval<-Eval(mu,alfa,beta) # el mejor valor
mejores<-matrix(c(mejor.eval,mu,alfa,beta,1),1,5)
k<-0 # indice de iteracion
actu<-0 # indice de actualizacion
set.seed(1) # Fijamos semilla para que siempre retorne el mismo resultado
plot(alfa,beta,xlim=c(0,2000),ylim=c(0,2000),type="n",xlab="alfa",ylab="beta",main = paste("Act:",actu,"Mejor alfa=",round(alfa,2),"Mejor beta=",round(beta,2)))
points(alfa,beta,cex=10,col="red",pch=3)

while (facred.acu>toler)
{
  k<-k+1
  # Genero nuevos valores aleatorios
  mu<-runif(1,mejor.mu-rango.mu*facred.acu,mejor.mu+rango.mu*facred.acu)
  alfa<-runif(1,mejor.alfa-rango.alfa*facred.acu,mejor.alfa+rango.alfa*facred.acu)
  beta<-runif(1,mejor.beta-rango.beta*facred.acu,mejor.beta+rango.beta*facred.acu)
  # Evaluacion de los nuevos valores
  valor<-Eval(mu,alfa,beta)
  if (valor<mejor.eval) # SI encuentro algo mejor -> Actualizacion
  {
    invisible(readline(prompt="Presione [enter] para seguir:"))
        actu<-actu+1
    # grafico
    plot(alfa,beta,xlim=c(0,2000),ylim=c(0,2000),type="n",xlab="alfa",ylab="beta",main =paste("Act:",actu,"Mejor alfa=",round(alfa,2),"Mejor beta=",round(beta,2)))
    points(alfa,beta,cex=10,col="red",pch=3)
    mejor.eval<-valor
    mejor.mu<-mu
    mejor.alfa<-alfa
    mejor.beta<-beta
    mejores<-rbind(mejores,c(mejor.eval,mu,alfa,beta,k))
  }
  else # SI NO encuentro algo mejor -> Reduzco rango de busqueda
  {
    facred.acu<-facred.acu*facred
  }
  points(alfa,beta,cex=0.1,col="blue") # puntos evalyados
}

NA

NA

NA

NA

NA

NA

NA

NA

NA

NA

NA

NA

NA

NA

NA

NA

NA

NA

NA

NA

NA

NA

NA

NA

NA

NA

NA

NA

NA

NA

NA

NA

NA

NA

NA

NA

NA

NA

NA

NA

c(mu,alfa,beta)

[1] 133400.6267    863.4835    590.5030

dim(mejores)

[1] 43  5

graficacion

plot((mejores[,1]),type="l",col="blue",xlab="Actualizacion",ylab="Valor de Perdida")
#lines(evol2[subconj,2],col="green")
abline(coe[1],0,lty=3)

title("Convergencia de la Perdida")

#
plot(mejores[,2],type="l",col="blue",xlab="Actualizacion",ylab="Valor de mu")
#lines(evol2[subconj,2],col="green")
abline(coe[1],0,lty=3)

title("Convergencia del Parametro mu")

#
#
plot(mejores[,3],type="l",col="blue",xlab="Actualizacion",ylab="Valor de alfa")
#lines(evol2[subconj,3],col="green")
abline(coe[2],0,lty=3)

title("Convergencia del Parametro alfa")

#
#
plot(mejores[,4],type="l",col="blue",xlab="Actualizacion",ylab="Valor de beta")
#lines(evol2[subconj,4],col="green")
abline(coe[3],0,lty=3)

title("Convergencia del Parametro beta")

#
#
plot(mejores[,5],type="l",col="blue",xlab="Actualizacion",ylab="Iteraciones")
#lines(evol2[subconj,4],col="green")
title("Evolucion de los Iteraciones")

#

¿ Cuanbueno es el ajuste ?

predichos<-predict(ajusM2)
mean(abs(Precio-predichos))

[1] 60757.25

pmaeM2<-mean(abs(Precio-predichos))/mean(Precio)
pmaeM2

[1] 0.1745913

Medida del error por Validacion Cruzada

n<-length(Precio)
predichos.oos<-rep(NA,n) # predichos out of sample
plot3d(ajusM2,size=15,col="blue")
invisible(readline(prompt="Presione [enter] para seguir:"))

for (i in 1:n)
{
  ajus.cv<-lm(Precio~Sup+Fondo,data=datos1d[-i,])
  predichos.oos[i]<-predict(ajus.cv,newdata=datos1d[i,])
  plot3d(ajus.cv,size=15,col="green",add=T,alpha=0.1)
}

# MAE
mean(abs(Precio-predichos.oos))

[1] 64086.7

# PMAE
pmaeM2.cv<-mean(abs(Precio-predichos.oos))/mean(Precio)
pmaeM2.cv

[1] 0.1841587

¿ Cuán certeras (creibles/estables/repetibles) son las relaciones halladas ? El Bootstrap

B<-1000 # cantidad de muestras bootstrap
mues<-rep(NA,B) # vector para guardar los mu estimados
alfas<-rep(NA,B) # vector para guardar los alfa estimados
betas<-rep(NA,B) # vector para guardar los beta estimados
set.seed(1)
for (b in 1:B)
{
  indices<-sample(1:68,68,replace = TRUE)
  ajus.boot<-lm(Precio~Sup+Fondo,data=datos1d[indices,])
  coe<-coef(ajus.boot)
  mues[b]<-coe[1]
  alfas[b]<-coe[2]
  betas[b]<-coe[3]
}
resul<-cbind(mues,alfas,betas)
head(resul)

         mues     alfas    betas
[1,] 150778.5  572.0181 910.9515
[2,] 144559.6  704.4977 764.7104
[3,] 142702.9  729.8691 606.0145
[4,] 172381.2  602.5592 679.1868
[5,] 115306.3  988.1902 660.7135
[6,] 109073.9 1023.9048 411.6738

Grafico de resultados

plot(alfas,betas)
# parametros estimados en el modelo inicial
segments(863,0,863,2000,col="green",lty = 1,lwd=3)

segments(0,591,2000,591,col="green",lty = 1,lwd=3)
# cuantiles bootstrap
segments(quantile(alfas,0.05),0,quantile(alfas,0.05),2000,col="red",lty = 3,lwd=3)

segments(quantile(alfas,0.5),0,quantile(alfas,0.5),2000,col="red",lty = 3,lwd=3)
segments(quantile(alfas,0.95),0,quantile(alfas,0.95),2000,col="red",lty = 3,lwd=3)

segments(0,quantile(betas,0.05),2000,quantile(betas,0.05),2000,col="red",lty = 3,lwd=3)
segments(0,quantile(betas,0.5),2000,quantile(betas,0.5),2000,col="red",lty = 3,lwd=3)

segments(0,quantile(betas,0.95),2000,quantile(betas,0.95),2000,col="red",lty = 3,lwd=3)
abline(0,1,col="blue")

title("Estimaciones Bootstrap de los Parametros")

mean(betas>alfas)

[1] 0.167

Modelo simple con Sup

ajusM1<-lm(Precio~Sup,data=datos1d)
coe<-coef(ajusM1)
coe

(Intercept)         Sup 
 109601.840    1254.226 

summary(ajusM1)

Call:
lm(formula = Precio ~ Sup, data = datos1d)

Residuals:
    Min      1Q  Median      3Q     Max 
-217225  -49083     603   75104  190740 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 109601.8    32593.3   3.363  0.00129 ** 
Sup           1254.2      160.6   7.808 5.83e-11 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 94100 on 66 degrees of freedom
Multiple R-squared:  0.4802,    Adjusted R-squared:  0.4723 
F-statistic: 60.97 on 1 and 66 DF,  p-value: 5.832e-11

plot(ajusM1)

plot(Sup,Precio)
abline(ajusM1)

title("Regresion Lineal Simple Precio Vs. Sup")

¿ Cuan bueno es el ajuste ?

predichos<-predict(ajusM1)
mean(abs(Precio-predichos))

[1] 76741.89

pmaeM1<-mean(abs(Precio-predichos))/mean(Precio)
pmaeM1

[1] 0.2205245

# cv
n<-length(Precio)
predichos.oos<-rep(NA,n) # predichos out of sample
for (i in 1:n)
{
  ajus.cv<-lm(Precio~Sup,data=datos1d[-i,])
  predichos.oos[i]<-predict(ajus.cv,newdata=datos1d[i,])
}
mean(abs(Precio-predichos.oos))

[1] 79186.78

pmaeM1.cv<-mean(abs(Precio-predichos.oos))/mean(Precio)
pmaeM1.cv

[1] 0.2275501

Modelo mas complejo

# agrego lat y lon
datos1e<-datos1a %>% filter(l3=="Boedo",property_type=="Casa")  %>% mutate(Precio=price,Sup=surface_covered,Fondo=surface_total-surface_covered) %>% dplyr::select(Sup,Fondo,Precio,lat,lon) %>%  filter(Fondo>=0) %>% na.omit()
datos1e

summary(datos1e)

      Sup            Fondo            Precio            lat              lon        
 Min.   : 35.0   Min.   :  0.00   Min.   : 18000   Min.   :-34.65   Min.   :-58.43  
 1st Qu.:132.0   1st Qu.:  7.00   1st Qu.:246500   1st Qu.:-34.64   1st Qu.:-58.42  
 Median :180.0   Median : 33.00   Median :297000   Median :-34.63   Median :-58.42  
 Mean   :185.7   Mean   : 82.56   Mean   :335166   Mean   :-34.63   Mean   :-58.41  
 3rd Qu.:251.0   3rd Qu.:156.50   3rd Qu.:450000   3rd Qu.:-34.62   3rd Qu.:-58.41  
 Max.   :332.0   Max.   :368.00   Max.   :590000   Max.   :-34.60   Max.   :-58.38  

ggpairs(datos1e)

# Creo las nuevas variables
Sup<-datos1e$Sup
Fondo<-datos1e$Fondo
Lon<-datos1e$lon
Lat<-datos1e$lat
Precio<-datos1e$Precio
ajusM4<-lm(Precio~Sup+Fondo+Lon+Lat)
coe<-coef(ajusM4)
coe

 (Intercept)          Sup        Fondo          Lon          Lat 
1.922485e+08 9.948234e+02 1.561275e+02 1.693947e+06 2.690522e+06 

summary(ajusM4)

Call:
lm(formula = Precio ~ Sup + Fondo + Lon + Lat)

Residuals:
    Min      1Q  Median      3Q     Max 
-202730  -37972      93   31036  217949 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 1.922e+08  5.647e+07   3.404  0.00113 ** 
Sup         9.948e+02  1.472e+02   6.758 4.35e-09 ***
Fondo       1.561e+02  1.606e+02   0.972  0.33453    
Lon         1.694e+06  1.372e+06   1.235  0.22122    
Lat         2.691e+06  1.712e+06   1.571  0.12087    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 78160 on 66 degrees of freedom
Multiple R-squared:  0.7103,    Adjusted R-squared:  0.6928 
F-statistic: 40.46 on 4 and 66 DF,  p-value: < 2.2e-16

plot(Lon,Lat)

predichos<-predict(ajusM4)
mean(abs(Precio-predichos))

[1] 53247.78

pmaeM4<-mean(abs(Precio-predichos))/mean(Precio)
pmaeM4

[1] 0.1588698

plot(Precio,predichos,xlab="Precios Observados",yla="Predichos por M4",main="Precios Predichos por M4 Vs. Observados")
abline(0,1,col="green")

Error por CV del modelo M4

n<-length(Precio)
predichos.oos<-rep(NA,n) # predichos out of sample
for (i in 1:n)
{
  ajus.cv<-lm(Precio~Sup+Fondo+lon+lat,data=datos1e[-i,])
  predichos.oos[i]<-predict(ajus.cv,newdata=datos1e[i,])
}
mean(abs(Precio-predichos.oos))

[1] 57397.51

pmaeM4.cv<-mean(abs(Precio-predichos.oos))/mean(Precio)
pmaeM4.cv

[1] 0.1712509

Vamos por mas …. Modelo mas mas complejo

agrego efectos cuadraticos a lat y lon

# agrego efectos no lineales
ajusM8<-lm(Precio~Sup+Fondo+poly(Lon,3)+poly(Lat,3))
coe<-coef(ajusM8)
coe

  (Intercept)           Sup         Fondo poly(Lon, 3)1 poly(Lon, 3)2 poly(Lon, 3)3 
  136293.1124      986.0497      190.8774   557146.4452   110893.3951   -42545.1995 
poly(Lat, 3)1 poly(Lat, 3)2 poly(Lat, 3)3 
  -76485.5627  -154553.0046   -78618.7351 

summary(ajusM8)

Call:
lm(formula = Precio ~ Sup + Fondo + poly(Lon, 3) + poly(Lat, 
    3))

Residuals:
    Min      1Q  Median      3Q     Max 
-216390  -34175     698   31947  219714 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)    136293.1    30714.0   4.437 3.80e-05 ***
Sup               986.0      154.2   6.394 2.37e-08 ***
Fondo             190.9      193.9   0.985    0.329    
poly(Lon, 3)1  557146.4  1591045.5   0.350    0.727    
poly(Lon, 3)2  110893.4   443020.8   0.250    0.803    
poly(Lon, 3)3  -42545.2   106206.2  -0.401    0.690    
poly(Lat, 3)1  -76485.6  1507667.2  -0.051    0.960    
poly(Lat, 3)2 -154553.0   761154.1  -0.203    0.840    
poly(Lat, 3)3  -78618.7   274829.6  -0.286    0.776    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 80210 on 62 degrees of freedom
Multiple R-squared:  0.7134,    Adjusted R-squared:  0.6764 
F-statistic: 19.29 on 8 and 62 DF,  p-value: 3.387e-14

plot(Lon,Lat)

predichos<-predict(ajusM8)
mean(abs(Precio-predichos))

[1] 52459.34

pmaeM8<-mean(abs(Precio-predichos))/mean(Precio)
pmaeM8

[1] 0.1565174

plot(Precio,predichos,xlab="Precios Observados",yla="Predichos por M8",main="Precios Predichos por M8 Vs. Observados")
abline(0,1,col="green")

Medida del error por Validacion Cruzada

n<-length(Precio)
predichos.oos<-rep(NA,n) # predichos out of sample
for (i in 1:n)
{
  ajus.cv<-lm(Precio~Sup+Fondo+poly(lon,3)+poly(lat,3),data=datos1e[-i,])
  predichos.oos[i]<-predict(ajus.cv,newdata=datos1e[i,])
}
mean(abs(Precio-predichos.oos))

[1] 61888.82

pmaeM8.cv<-mean(abs(Precio-predichos.oos))/mean(Precio)
pmaeM8.cv

[1] 0.1846511

Comparacion de Modelos

PMAES<-matrix(c(pmaeM1,pmaeM1.cv,pmaeM2,pmaeM2.cv,pmaeM4,pmaeM4.cv,pmaeM8,pmaeM8.cv),4,2,byrow = T)
PMAES

          [,1]      [,2]
[1,] 0.2205245 0.2275501
[2,] 0.1745913 0.1841587
[3,] 0.1588698 0.1712509
[4,] 0.1565174 0.1846511

matplot(PMAES,type="b",pch=16,col=2:3,lty=1,lwd=3,, xaxt='n',main="PMAEs Ingenuos y PMAEs CV para M!, M2, M4 y M8",ylim=c(0.13,0.25),xlab="Modelo",ylab="PMAE")
text(1,pmaeM1.cv+0.01,"M1")

text(2,pmaeM2.cv+0.01,"M2")
text(3,pmaeM4.cv+0.01,"M4")

text(4,pmaeM8.cv+0.01,"M8")

Trade-off Sesgo-Varianza

N<-68
Sup0<-250
Relacion<-function(x){(log((x-70)/520)-log(10/520))*100000+100000}
set.seed(1)
Sup<-runif(N)*520+80
Error<-rnorm(N)*50000
Precio<-Relacion(Sup)+Error
plot(Sup,Precio,col="blue",main="Precio Vs. Superficie")
curve(Relacion,from=80,to=600,add=TRUE,col="green",lwd=2)

Precio0<-Relacion(Sup0)
segments(Sup0,0,Sup0,Precio0,lty=3,col="green",lwd=2)
segments(0,Precio0,Sup0,Precio0,lty=3,col="green",lwd=2)

Estimaciones simples

set.seed(1)
Sup<-runif(N)*520+80
Error<-rnorm(N)*50000
Precio<-Relacion(Sup)+Error
cant<-50 # cantidad de estimaciones
estimaciones<-rep(NA,cant)
grilla<-data.frame(Sup=seq(80,600,length.out = 100))
# grafico fijo
plot(Sup,Precio,col="blue",main="Precio Vs. Superficie con Ajustes de Modelos Simples")

# bucle
for (i in 1:cant)
{
Sup<-runif(N)*520+80
Error<-rnorm(N)*50000
Precio<-Relacion(Sup)+Error
ajus.1<-lm(Precio~poly(Sup,1))
pred.1<-predict(ajus.1,newdata=data.frame(Sup=Sup0))
estimaciones[i]<-pred.1
curva.1<-predict(ajus.1,newdata = grilla)
lines(as.numeric(grilla$Sup),curva.1,lwd=0.5,col="red")
}

segments(0,mean(estimaciones),Sup0,mean(estimaciones),lty=3,lwd=2,col="black")
segments(0,mean(estimaciones)-2*sd(estimaciones),Sup0,mean(estimaciones)-2*sd(estimaciones),lty=3,lwd=2,col="black")

segments(0,mean(estimaciones)+2*sd(estimaciones),Sup0,mean(estimaciones)+2*sd(estimaciones),lty=3,lwd=2,col="black")
#
curve(Relacion,from=80,to=600,add=TRUE,col="green",lwd=3)

segments(Sup0,0,Sup0,Precio0,lty=3,lwd=2,col="green")
segments(0,Precio0,Sup0,Precio0,lty=3,lwd=2,col="green")

estimaciones.simples<-estimaciones

Estimaciones complejas

set.seed(1)
Sup<-runif(N)*520+80
Error<-rnorm(N)*50000
Precio<-Relacion(Sup)+Error
cant<-50 # cantidad de estimaciones
estimaciones<-rep(NA,cant)
grilla<-data.frame(Sup=seq(80,600,length.out = 100))
# grafico fijo
plot(Sup,Precio,col="blue",main="Precio Vs. Superficie con Ajustes de Modelos Complejos")

# bucle
for (i in 1:cant)
{
Sup<-runif(N)*520+80
Error<-rnorm(N)*50000
Precio<-Relacion(Sup)+Error
ajus.1<-lm(Precio~poly(Sup,12))
pred.1<-predict(ajus.1,newdata=data.frame(Sup=Sup0))
estimaciones[i]<-pred.1
curva.1<-predict(ajus.1,newdata = grilla)
lines(as.numeric(grilla$Sup),curva.1,lwd=0.5,col="red")
}

segments(0,mean(estimaciones),Sup0,mean(estimaciones),lty=3,lwd=2,col="black")
segments(0,mean(estimaciones)-2*sd(estimaciones),Sup0,mean(estimaciones)-2*sd(estimaciones),lty=3,lwd=2,col="black")

segments(0,mean(estimaciones)+2*sd(estimaciones),Sup0,mean(estimaciones)+2*sd(estimaciones),lty=3,lwd=2,col="black")#
curve(Relacion,from=80,to=600,add=TRUE,col="green",lwd=3)

segments(Sup0,0,Sup0,Precio0,lty=3,lwd=2,col="green")
segments(0,Precio0,Sup0,Precio0,lty=3,lwd=2,col="green")

estimaciones.complejas<-estimaciones

Comparacion

  boxplot(list(Simple=estimaciones.simples,Complejo=estimaciones.complejas),ylab="Precio",main="Distribucion de las Predicciones",col=c("cyan","magenta"))
abline(Precio0,0,col="green",lty=3,lwd=2)

Comparacion modelo M1 con M2

ajus.Fondo<-lm(Fondo~Sup,data=datos1e)
coe.Fondo<-coef(ajus.Fondo)
coe.M1<-coef(ajusM1)
coe.M2<-coef(ajusM2)
# pendiente
coe.M2[2]+coe.M2[3]*coe.Fondo[2]

     Sup 
1254.545 

# intercept
coe.M2[1]+coe.M2[3]*coe.Fondo[1]

(Intercept) 
   109533.2 

# el ajuste M1
coe.M1

(Intercept)         Sup 
 109601.840    1254.226