load("/Users/utente/Documents/prova r pkg/lucidi del corso/Kuznets.RData")
plot(provku$rmed,provku$gindex,cex=0.5,pch=19,xlim=c(13,28),ylim=c(0.36,0.48),xlab="Per Capita Income",ylab="Gindex",main="Provincial and Regional Per Capita Income and G-Index") # su questa relazione c'è una teoria - da Nobel - che dovrebbe essere parabolica, col cazzo !
points(regku$rmed,regku$gindex,cex=1,pch=19,col=2)
text(regku$rmed,regku$gindex,regku$cod,pos = 4,cex=1.2,col=2)
text(provku$rmed,provku$gindex,provku$cod,pos = 4)
abline(v=itku$rmed,lw=3)
abline(h=itku$gindex,lw=3)
text(16.5,0.48,"Per Capita Income in Italy",pos=4)
text(26,0.43,"G-Index in Italy",pos=4)

library(MASS)
plot(provku$rmed,provku$gindex,cex=0.5,pch=19,xlim=c(13,28),ylim=c(0.36,0.48),xlab="Per Capita Income",ylab="Gindex",main="Provincial and Regional Per Capita Income and G-Index")
reglm   <-lm(gindex~rmed, data=provku)
regrml1 <-rlm(gindex~rmed, data=provku, psi=psi.huber) 
regrml2 <-rlm(gindex~rmed, data=provku, psi=psi.bisquare) 
reglts  <-lqs(gindex~rmed, data=provku, method="lts")
abline(reglm)
abline(regrml1,lwd=2)
abline(regrml2,lwd=3)
abline(reglts,lwd=4)
confronto<-data.frame(coef(reglm),coef(regrml1),coef(regrml2),coef(reglts))
options(digits=3)
confronto

library(car)
outlierTest(reglm)
lev<-hat(model.matrix(reglm))
lev<-hatvalues(reglm)
n<-length(lev)
p<-sum(lev)
plot(lev, main="Punti di leva",cex=1,pch=19)
abline(h=2*p/n)

cook <- cooks.distance(reglm)
infl <- influence.measures(reglm)

load("/Users/utente/Documents/prova r pkg/lucidi del corso/lifexp.RData")
pairs(~Life.expectancy+Alcohol+percentage.expenditure+BMI+Polio+Total.expenditure+Diphtheria+log(GDP)+Income.composition.of.resources+Schooling,data = ledat[ledat$Year==2010,],cex=0.5,pch=19,main="Anno 2010")
library(MASS)
regfull <- lm(Life.expectancy ~ Year+Alcohol+percentage.expenditure+BMI+Polio+Total.expenditure+Diphtheria+log(GDP)+Income.composition.of.resources+Schooling, data = ledat)
ledat2010 <- ledat[ledat$Year==2010,]
ledat2010 <- ledat2010[complete.cases(ledat2010),]
regfull <- lm(Life.expectancy ~ Alcohol+percentage.expenditure+BMI+Polio+Total.expenditure+Diphtheria+log(GDP)+Income.composition.of.resources+Schooling, data = ledat2010)
regst   <- stepAIC(regfull, direction = "both",na.rm=T)
summary(regst)

swiss <- datasets::swiss
x <- model.matrix(Fertility~., swiss)[,-1]
y <- swiss$Fertility
lambda <- 10^seq(10, -2, length = 100)
swisslm <- lm(Fertility~., data = swiss)
coef(swisslm)
library(glmnet)
ridge.mod <- glmnet(x, y, alpha = 0, lambda = lambda)
predict(ridge.mod, s = 0, type = 'coefficients')[1:6,]
set.seed(489)
train = sample(1:nrow(x), nrow(x)/2)
test = (-train)
ytest = y[test]
swisslm <- lm(Fertility~., data = swiss, subset = train)
ridge.mod <- glmnet(x[train,], y[train], alpha = 0, lambda = lambda)
cv.out <- cv.glmnet(x[train,], y[train], alpha = 0)
bestlam <- cv.out$lambda.min
ridge.pred <- predict(ridge.mod, s = bestlam, newx = x[test,])
s.pred <- predict(swisslm, newdata = swiss[test,])
mean((s.pred-ytest)^2)
mean((ridge.pred-ytest)^2)
out = glmnet(x[train,],y[train],alpha = 0)
predict(ridge.mod, type = "coefficients", s = bestlam)[1:6,]
lasso.mod <- glmnet(x[train,], y[train], alpha = 1, lambda = lambda)
lasso.pred <- predict(lasso.mod, s = bestlam, newx = x[test,])
mean((lasso.pred-ytest)^2)
lasso.coef  <- predict(lasso.mod, type = 'coefficients', s = bestlam)[1:6,]

set.seed(999)
cv.ridge <- cv.glmnet(x, y, family='binomial', alpha=0, parallel=TRUE, standardize=TRUE, type.measure='auc')
plot(cv.ridge)
cv.ridge$lambda.min
cv.ridge$lambda.1se
coef(cv.ridge, s=cv.ridge$lambda.min)
set.seed(999)
cv.lasso <- cv.glmnet(x, y, family='binomial', alpha=1, parallel=TRUE, standardize=TRUE, type.measure='auc')
plot(cv.lasso)
plot(cv.lasso$glmnet.fit, xvar="lambda", label=TRUE)
cv.lasso$lambda.min
cv.lasso$lambda.1se
coef(cv.lasso, s=cv.lasso$lambda.min)

load("/Users/utente/Documents/prova r pkg/lucidi del corso/worldhappiness.RData")
plot(whdat$Economy..GDP.per.Capita.,whdat$Happiness.Score,cex=0.5,pch=19,main = "Felicità e Reddito",ylab="Indice di Felicità",xlab="GDP pro capite")
cor(whdat$Economy..GDP.per.Capita.,whdat$Happiness.Score)
hist(whdat$Happiness.Score)
reg <- lm(whdat$Happiness.Score ~ whdat$Economy..GDP.per.Capita.)
library(quantreg)
taus <- c(.05,.1,.25,.5,.75,.9,.95)
taus <- seq(0.05,0.95,0.1)
regq <- rq(whdat$Happiness.Score ~ whdat$Economy..GDP.per.Capita.,tau = taus)
plot(regq)
plot(whdat$Happiness.Score ~ whdat$Economy..GDP.per.Capita., data = whdat, cex=0.5,pch = 19, main = "Felicità e Reddito",ylab="Indice di Felicità",xlab="GDP pro capite")
for (j in 1:ncol(regq$coefficients)) {
  abline(coef(regq)[, j])
}
hist(ledat2010$Life.expectancy)
regq <- rq(ledat2010$Life.expectancy ~ ledat2010$Schooling,tau = taus)
plot(regq)
plot(ledat2010$Life.expectancy ~ ledat2010$Schooling, cex=0.5,pch = 19, main = "Salute e Scuola",xlab="Scolarizzazione",ylab="Aspettativa di Vita")
for (j in 1:ncol(regq$coefficients)) {
  abline(coef(regq)[, j])
}

load("/Users/utente/Documents/prova r pkg/lucidi del corso/penis.RData")
plot(penis$volume_flaccid,penis$volume_erect,cex=1,pch=19,main = "Penis Data",xlab="Volume a riposo",ylab="Volume eretto",col=penis$Region,cex.lab=2)
regns <- lm(penis$volume_erect ~ penis$volume_flaccid * penis$Region)
set.seed(2807)
library("flexmix")
Model <- FLXMRglm(~ volume_flaccid)
regmis <- stepFlexmix(volume_erect ~ 1, model = Model, nrep = 3, k = 5, data = penis, concomitant = FLXPmultinom(formula = ~1))
regmis <- relabel(regmis, "model", "volume_flaccid")
summary(refit(regmis))
library(gstat)
data(meuse.riv)
data(meuse.area)
par(mar=c(1,1,1,1),mfrow=c(1,2))
plot(meuse.area, type = "l", asp = 1,axes=F)
polygon(meuse.riv,col=5)
points(meuse.all[,2:3],cex=1,pch=19)
par(mar=c(4.5,4.5,1,1))
plot(meuse.all$lead,meuse.all$cadmium,cex=1,pch=19,main="",xlab="Piombo",ylab="Cadmio",cex.lab=2)
library(mixreg)
mixr <- mixreg(meuse.all$lead,meuse.all$cadmium,ncomp=2)
cvmr <- covmix(mixr,meuse.all$lead,meuse.all$cadmium)
cbdr <- cband(mixr,cvmr,meuse.all$lead,meuse.all$cadmium)
plot(cbdr)
load("/Users/utente/Documents/Siena 2016/mandalaz.RData")
par(mar=c(1,1,1,1),mfrow=c(1,1))
plot(zberg$x_terr,zberg$y_terr,cex=1,pch=19,axes=F)
par(mar=c(4.5,4.5,1,1),mfrow=c(1,2))
plot(zberg$stem,zberg$basal,cex=1,pch=19,main="",xlab="Gambo",ylab="Basale",cex.lab=2)
mixr <- mixreg(zberg$stem,zberg$basal,ncomp=2)
cvmr <- covmix(mixr,zberg$stem,zberg$basal)
cbdr <- cband(mixr,cvmr,zberg$stem,zberg$basal)
plot(cbdr)

load("/Users/utente/Documents/prova r pkg/lucidi del corso/HumanFreedomIndex.RData")
datisel <- hfi[,c(1,3,4,8,22,26,32,44,52,61,62,70,81,86,99,118,119,121)]
datis <- datisel[datisel$year==2010,]

library(Ecdat)
reg = lm(nbaffairs~sex+age+child+ym+religious+education+occupation+rate,data=Fair)
x=model.matrix(nbaffairs~sex+age+child+ym+religious+education+occupation+rate,Fair)[,-1]
y=Fair$nbaffairs
library(glmnet)
grid=10^seq(10,-2,length=100)
ridge.mod=glmnet(x,y,alpha=0,lambda=grid)
dim(coef(ridge.mod))
ridge.mod$lambda[50]
coef(ridge.mod)[,50]
sqrt(sum(coef(ridge.mod)[-1,50]^2))
ridge.mod$lambda[60]
coef(ridge.mod)[,60]
sqrt(sum(coef(ridge.mod)[-1,60]^2))
predict(ridge.mod,s=50,type="coefficients")[1:9,]
set.seed(1)
train=sample(1:nrow(x), nrow(x)/2)
test=(-train)
y.test=y[test]
ridge.mod=glmnet(x[train,],y[train],alpha=0,lambda=grid, thresh=1e-12)
ridge.pred=predict(ridge.mod,s=4,newx=x[test,])
mean((ridge.pred-y.test)^2)
mean((mean(y[train])-y.test)^2)
ridge.pred=predict(ridge.mod,s=1e10,newx=x[test,])
mean((ridge.pred-y.test)^2)
ridge.pred=predict(ridge.mod,s=0,newx=x[test,],exact=T,x=x[train,],y=y[train])
mean((ridge.pred-y.test)^2)
lm(y~x, subset=train)
predict(ridge.mod,s=0,exact=T,type="coefficients",x=x[train,],y=y[train])[1:9,]
set.seed(1)
cv.out=cv.glmnet(x[train,],y[train],alpha=0)
bestlam=cv.out$lambda.min
bestlam
ridge.pred=predict(ridge.mod,s=bestlam,newx=x[test,])
mean((ridge.pred-y.test)^2)
out=glmnet(x,y,alpha=0)
predict(out,type="coefficients",s=bestlam)[1:9,]
par(mfrow=c(1,2))
plot(cv.out, main="Ridge")

lasso.mod=glmnet(x[train,],y[train],alpha=1,lambda=grid)
#plot(lasso.mod)
set.seed(1)
cv.out=cv.glmnet(x[train,],y[train],alpha=1)
plot(cv.out,main="Lasso")
bestlam=cv.out$lambda.min
lasso.pred=predict(lasso.mod,s=bestlam,newx=x[test,])
mean((lasso.pred-y.test)^2)
out=glmnet(x,y,alpha=1,lambda=grid)
lasso.coef=predict(out,type="coefficients",s=bestlam)[1:9,]
lasso.coef
lasso.coef[lasso.coef!=0]

library(Ecdat)
lreg <- lm(price~carat,data=Diamond)
fit=lm(price~poly(carat,3),data=Diamond)
coef(summary(fit))
fit2=lm(price~poly(carat,3,raw=T),data=Diamond)
coef(summary(fit2))
fit2a=lm(price~carat+I(carat^2)+I(carat^3),data=Diamond)
coef(fit2a)
fit2b=lm(price~cbind(carat,carat^2,carat^3),data=Diamond)
caratlims=range(Diamond$carat)
carat.grid=seq(from=caratlims[1],to=caratlims[2],by=0.01)
preds=predict(fit,newdata=list(carat=carat.grid),se=TRUE)
se.bands=cbind(preds$fit+2*preds$se.fit,preds$fit-2*preds$se.fit)
par(mfrow=c(1,1),mar=c(4.5,4.5,1,1),oma=c(0,0,4,0))
plot(Diamond$carat,Diamond$price,xlim=caratlims,cex=.5,pch=19,col="darkgrey",xlab="Carati",ylab="Prezzo",main="Degree-3 Polynomial",cex.lab=2)
lines(carat.grid,preds$fit,lwd=2,col="blue")
matlines(carat.grid,se.bands,lwd=1,col="blue",lty=3)
abline(lreg,lwd=2,col=2)
preds2=predict(fit2,newdata=list(carat=carat.grid),se=TRUE)
max(abs(preds$fit-preds2$fit))
fit.1=lm(price~carat,data=Diamond)
fit.2=lm(price~poly(carat,2),data=Diamond)
fit.3=lm(price~poly(carat,3),data=Diamond)
fit.4=lm(price~poly(carat,4),data=Diamond)
fit.5=lm(price~poly(carat,5),data=Diamond)
anova(fit.1,fit.2,fit.3,fit.4,fit.5)
Diamond <- Diamond[order(Diamond$carat),]
loess25 <- loess(price ~ carat, data=Diamond, span=0.25)
loess50 <- loess(price ~ carat, data=Diamond, span=0.50)
loess75 <- loess(price ~ carat, data=Diamond, span=0.75)
smoothed25 <- predict(loess25) 
smoothed50 <- predict(loess50) 
smoothed75 <- predict(loess75) 
plot(Diamond$carat, Diamond$price, main="Loess Smoothing and Prediction", xlab="Carati", ylab="Prezzo",xlim=caratlims,cex=.5,pch=19,col="darkgrey",cex.lab=2)
lines(smoothed25, x=Diamond$carat, col="red",lwd=2)
lines(smoothed50, x=Diamond$carat, col="blue",lwd=2)
lines(smoothed75, x=Diamond$carat, col="black",lwd=2)
abline(lreg,lwd=2,col="green")

load("/Users/utente/Documents/prova r pkg/lucidi del corso/esoplan.RData")
esop <- esopianeti[!is.na(esopianeti$MASS) & !is.na(esopianeti$RADIUS),]
esop <- esop[order(esop$MASS),]
lmas <- log(esop$MASS)
lrag <- log(esop$RADIUS)
loess25 <- loess(lrag ~ lmas, span=0.25)
loess50 <- loess(lrag ~ lmas, span=0.50)
loess75 <- loess(lrag ~ lmas, span=0.75)
smoothed25 <- predict(loess25) 
smoothed50 <- predict(loess50) 
smoothed75 <- predict(loess75) 
plot(log(esop$MASS),log(esop$RADIUS),cex=0.5,pch=19,xlab="Massa", ylab="Raggio",cex.lab=2)
lines(smoothed25, x=lmas, col="red",lwd=2)
lines(smoothed50, x=lmas, col="blue",lwd=2)
lines(smoothed75, x=lmas, col="black",lwd=2)
load("/Users/utente/Documents/prova r pkg/lucidi del corso/fifa19.RData")
fifa <- fifa19[order(fifa19$Overall),]
loess25 <- loess(Valore ~ Overall, data=fifa, span=0.25)
loess50 <- loess(Valore ~ Overall, data=fifa, span=0.50)
loess75 <- loess(Valore ~ Overall, data=fifa, span=0.75)
smoothed25 <- predict(loess25) 
smoothed50 <- predict(loess50) 
smoothed75 <- predict(loess75) 
plot(fifa$Overall,fifa$Valore,pch=1,cex=fifa$International.Reputation/3,xlab="Capacità", ylab="Valore",cex.lab=2)
lines(smoothed25, x=fifa$Overall, col="red",lwd=2)
lines(smoothed50, x=fifa$Overall, col="blue",lwd=2)
lines(smoothed75, x=fifa$Overall, col="black",lwd=2)
setwd("~/Documents/prova r pkg/lucidi del corso/Life Expectancy Data")
ledat <- read.csv("Life Expectancy Data.csv")[,c(4,17)]
ledat <- ledat[complete.cases(ledat),]
ledat <- ledat[order(ledat$GDP),]
loess25 <- loess(Life.expectancy ~ GDP, data=ledat, span=0.25)
loess50 <- loess(Life.expectancy ~ GDP, data=ledat, span=0.50)
loess75 <- loess(Life.expectancy ~ GDP, data=ledat, span=0.75)
smoothed25 <- predict(loess25) 
smoothed50 <- predict(loess50) 
smoothed75 <- predict(loess75) 
plot(ledat$GDP,ledat$Life.expectancy,pch=19,cex=0.5,xlab="GDP", ylab="Aspettativa di vita",cex.lab=2)
lines(smoothed25, x=ledat$GDP, col="red",lwd=2)
lines(smoothed50, x=ledat$GDP, col="blue",lwd=2)
lines(smoothed75, x=ledat$GDP, col="black",lwd=2)

library(splines)
fit=lm(price~bs(carat,knots=c(0.35,0.60,0.85,1,00)),data=Diamond)
pred=predict(fit,newdata=list(carat=carat.grid),se=T)
par(mfrow=c(1,2),mar=c(4.5,4.5,1,1))
plot(Diamond$carat,Diamond$price,pch=19,cex=0.5,col="gray",xlab="Carati",ylab="Prezzo",main="Spline Regression",cex.lab=2)
lines(carat.grid,pred$fit,lwd=2)
lines(carat.grid,pred$fit+2*pred$se,lty="dashed")
lines(carat.grid,pred$fit-2*pred$se,lty="dashed")
abline(v=0.35,col=2,lty=2)
abline(v=0.60,col=2,lty=2)
abline(v=0.85,col=2,lty=2)
abline(v=1.00,col=2,lty=2)
dim(bs(Diamond$carat,knots=c(0.35,0.60,0.85,1,00)))
dim(bs(Diamond$carat,df=6))
attr(bs(Diamond$carat,df=6),"knots")
fit2=lm(price~ns(carat,df=4),data=Diamond)
pred2=predict(fit2,newdata=list(carat=carat.grid),se=T)
plot(Diamond$carat,Diamond$price,xlim=caratlims,pch=19,cex=.5,col="darkgrey",,xlab="Carati",ylab="Prezzo",main="Smoothing Spline",cex.lab=2)
lines(carat.grid, pred2$fit,col="red",lwd=2)
fit=smooth.spline(Diamond$carat,Diamond$price,df=16)
fit2=smooth.spline(Diamond$carat,Diamond$price,cv=TRUE)
fit2$df
lines(fit,col="green",lwd=2)
lines(fit2,col="blue",lwd=2)
legend("topright",legend=c("16 DF","6.8 DF"),col=c("green","blue"),lty=1,lwd=2,cex=.8)

fit2=lm(lrag~ns(lmas,df=4),data=esop)
lmaslims=range(lmas)
lmas.grid=seq(from=lmaslims[1],to=lmaslims[2],by=0.01)
pred2 <- predict(fit2,newdata=list(lmas=lmas.grid),se=T)
par(mfrow=c(1,1),mar=c(4.5,4.5,1,1))
plot(log(esop$MASS),log(esop$RADIUS),cex=0.5,pch=19,xlab="Massa", ylab="Raggio",cex.lab=2)
lines(pred2$fit, x=lmas.grid, col="red",lwd=2)
lines(pred2$fit+1.96*pred2$se.fit, x=lmas.grid, col="blue",lwd=2,lty=2)
lines(pred2$fit-1.96*pred2$se.fit, x=lmas.grid, col="blue",lwd=2,lty=2)

fit2=lm(Valore~ns(Overall,df=4),data=fifa)
Overalllims=range(fifa$Overall)
Overall.grid=seq(from=Overalllims[1],to=Overalllims[2],by=0.1)
pred2 <- predict(fit2,newdata=list(Overall=Overall.grid),se=T)
par(mfrow=c(1,1),mar=c(4.5,4.5,1,1))
plot(fifa$Overall,fifa$Valore,pch=1,cex=fifa$International.Reputation/3,xlab="Capacità", ylab="Valore",cex.lab=2)
lines(pred2$fit, x=Overall.grid, col="red",lwd=4)

fit2=lm(Life.expectancy~ns(GDP,df=4),data=ledat)
GDPlims=range(ledat$GDP)
GDP.grid=seq(from=GDPlims[1],to=GDPlims[2],by=100)
pred2 <- predict(fit2,newdata=list(GDP=GDP.grid),se=T)
par(mfrow=c(1,1),mar=c(4.5,4.5,1,1))
plot(ledat$GDP,ledat$Life.expectancy,pch=19,cex=0.5,xlab="GDP", ylab="Aspettativa di vita",cex.lab=2)
lines(pred2$fit, x=GDP.grid, col="red",lwd=4)
lines(pred2$fit+1.96*pred2$se.fit, x=GDP.grid, col="blue",lwd=2,lty=2)
lines(pred2$fit-1.96*pred2$se.fit, x=GDP.grid, col="blue",lwd=2,lty=2)

gamlm <- lm(cons~ns(income,3)+ns(price,3)+ns(temp,3),data=Icecream)
library(gam)
gamice <- gam(cons~ns(income,3)+ns(price,3)+ns(temp,3),data=Icecream)
par(mfrow=c(1,3))
plot(gamice, se=TRUE,col="blue")
summary(gamice)

gamlm <- lm(Valore~ns(Overall,3)+ns(Age,3)+ns(Clausola,3)+International.Reputation,data=fifa)
library(gam)
gamfifa <- gam(Valore~ns(Overall,3)+ns(Age,3)+ns(Clausola,3)+International.Reputation,data=fifa)
par(mfrow=c(1,4))
plot(gamfifa, se=TRUE,col="blue")
summary(gamfifa)

load("/Users/utente/Documents/prova r pkg/lucidi del corso/movies.RData")
regfilm <- lm(gross ~ ns(budget,3) + ns(score,3) + ns(year,3) + ns(votes,3),data=revfilm)
library(gam)
gamfilm <- gam(gross ~ ns(budget,3) + ns(score,3) + ns(year,3) + ns(votes,3),data=revfilm)
par(mfrow=c(1,4))
plot(gamfilm, se=TRUE,col="blue")
summary(gamfilm)

library(AER)
data(OECDGrowth)
solow_lm <- lm(log(gdp85/gdp60) ~ log(gdp60)+log(invest)+log(popgrowth+.05)+log(school)+log(randd),data=OECDGrowth)
summary(solow_lm)
par(mfrow=c(1,1),mar=c(4.5,4.5,1,1))
plot(log(OECDGrowth$gdp60),log(OECDGrowth$gdp85/OECDGrowth$gdp60),cex=1,pch=19,xlab="log(GDP60)",ylab="log(GDP85/GDP60)",cex.lab=2)
text(log(OECDGrowth$gdp60),log(OECDGrowth$gdp85/OECDGrowth$gdp60),rownames(OECDGrowth),pos=3,cex=1.5)
pp <- (1/OECDGrowth$gdp60)/sum(1/OECDGrowth$gdp60)
set.seed(4000)
OECDGrowth_miss <- OECDGrowth
OECDGrowth_miss[sample(1:nrow(OECDGrowth),3,replace=F,prob=pp),3] <- NA
OECDGrowth_miss[sample(1:nrow(OECDGrowth),3,replace=F,prob=pp),4] <- NA
OECDGrowth_miss[sample(1:nrow(OECDGrowth),3,replace=F,prob=pp),5] <- NA
OECDGrowth_miss[sample(1:nrow(OECDGrowth),3,replace=F,prob=pp),6] <- NA
solow_lm_miss <- lm(log(gdp85/gdp60) ~ log(gdp60)+log(invest)+log(popgrowth+.05)+log(school)+log(randd),data=OECDGrowth_miss)
summary(solow_lm_miss)
library(miceadds)
md.pattern(OECDGrowth_miss, plot = T)
OECDGrowth_imp <- mice(OECDGrowth_miss, m = 15, seed = 4000)
lmimp <- with(OECDGrowth_imp, lm(log(gdp85/gdp60) ~ log(gdp60)+log(invest)+log(popgrowth+.05)+log(school)+log(randd)))
clmimp <- summary(pool(lmimp))
library(mvdalab)
# OECDGrowth_miss <- introNAs(OECDGrowth, percent = 25)
OECDGrowth_EM <- imputeEM(OECDGrowth_miss,impute.ncomps = 1)
solow_lm_EM <- lm(log(gdp85/gdp60) ~ log(gdp60)+log(invest)+log(popgrowth+.05)+log(school)+log(randd),data=as.data.frame(OECDGrowth_EM$Imputed.DataFrames))
OECDGrowth_BS <- imputeBasic(OECDGrowth_miss)
solow_lm_BS <- lm(log(gdp85/gdp60) ~ log(gdp60)+log(invest)+log(popgrowth+.05)+log(school)+log(randd),data=OECDGrowth_BS$Imputed.DataFrame)
compare <- cbind(summary(solow_lm)$coefficients[,1],summary(solow_lm_miss)$coefficients[,1],summary(solow_lm_BS)$coefficients[,1],summary(solow_lm_EM)$coefficients[,1],clmimp$estimate)
compstd <- cbind(summary(solow_lm)$coefficients[,2],summary(solow_lm_miss)$coefficients[,2],summary(solow_lm_BS)$coefficients[,2],summary(solow_lm_EM)$coefficients[,2],clmimp$std.error)

load("/Users/utente/Documents/prova r pkg/lucidi del corso/computer.RData")
hdati[1,2] <- 6
hdati[274,3] <- 30
data <- as.Date(ISOdate(hdati$anno,hdati$mm,hdati$gg))
lncap <- log(hdati$cap)
plot(data,lncap,cex=0.5,pch=19,xlab="Anni",ylab="Log(Capacità Hard-Disk)",main="Funzione sigmoidale per sviluppo di Moore",cex.lab=1.5)
moore <- nls(lncap~phi0+phi1/(1+exp(-(phi2+phi3*as.numeric(data)))),start=list(phi0=-5,phi1=15,phi2=-10,phi3=0.001),trace=TRUE)
pred <- predict(moore,se.fit = T)
lines(as.numeric(data),pred,lw=4,col="red")

plot(mdati$anno,mdati$transistors,cex=1,pch=19,xlab="Anni",ylab="Transistor",main="Prima legge di Moore",cex.lab=1.5)
moore <- nls(transistors~exp(b0)*exp(b1*anno),data=mdati,start=list(b0=1,b1=0.2),trace=TRUE,nls.control(maxiter=1000))
summary(moore)
summary(rr <- lm(log(transistors)~anno,data=mdati))
pred <- predict(moore,se.fit = T)
lines(mdati$anno,pred,lw=4,col="red")
1/(coef(moore)[2]/log(2)); 1/(coef(rr)[2]/log(2))

library(micEconCES)
data( germanFarms, package = "micEcon" )
germanFarms$qOutput <- germanFarms$vOutput / germanFarms$pOutput
germanFarms$qVarInput <- germanFarms$vVarInput / germanFarms$pVarInput
cesLandLabor <- cesEst( "qOutput", c( "land", "qLabor" ), germanFarms )
cesLandLaborNm <- cesEst( "qOutput", c( "land", "qLabor" ), germanFarms,method = "NM" )
data("GermanIndustry")
cesInd <- cesEst( "Y", c("K","A","E"), GermanIndustry,method = "NM")
summary(cesInd)