problemi di allocazione di memoria
Italian is not spoken on this list. And you were asked not to post in html. You did and your code is totally messed up and thus completely unreadable. In addition it is too much. Berend
On 17-06-2013, at 09:22, laura bettella <cyri at live.it> wrote:
Buongiorno,
io sto usando R per fare un lavoro di simulazione per la mia tesi. Il lavoro consiste nel simulare 100 modelli da replicare 300 volte, fare le previsione, e combinarle usando determinati pesi.Il problema purtroppo ? che il programma di ferma a 30 modelli perch? dice non non avere abbastanza spazio nella memoria. Pensavo di aver ovviato a questo problema andando in: risorse del computer -> C -> propriet? -> gestione quote e mettendo che non ci sia un limite alle quote. Lo spazio disponibile in C ? di 417 GB.
Di seguito vi metto il mio codice (so che pu? sembrare lungo):
# liberiamo l'ambiente di lavoro e carichiamo i pacchetti utilirm(list=objects())library(VGAM) # serve per generare dati dall'esponenziale doppia (Laplace)
# percorsi e parametri del codicebasepath = "C:/Documents and Settings/Laura/Desktop/top secret -)/lavoroLaura/lavoroLaura/"outputpath = paste0(basepath, "outputbreak/")nrep = 100 # numero replicazioni (numero modelli AR2 da cui simuleremo)#nrep = 1 # per provansim = 300 # numero simulazioni per ogni modello#nsim = 3 # per provaseme = 123 # seme simulazione dati (per replicare i risultati ottenuti)set.seed(seme)
# varianza sigma2 = 3
# funzione phi_s(.) di H-AFTERphis = function(x,s) { if (x>=-1 & x<=s) { x^2 } else { if (x>s) { 2*s*x-s^2 } else { -2*x-1 } }}
# s e lambda: parametri (fissi) di H-AFTER e L1-AFTERs = 1lambda = 1
# dataframe che conterr? gli nrep vettori di parametri degli AR3 da cui simuleremods_coeff_ar25 = data.frame(theta1=double(nrep*50), theta2=double(nrep*50), theta3=double(nrep*50), theta4=double(nrep*50), theta5=double(nrep*50))
# codice per tenere solo le combinazioni di parametri che danno luogo a modelli stazionari# (devono essere stazionari sia i modelli ar2, che gli ar5)for (i in 1:(nrep*50)) { coeff_ar12 = runif(2,-1,1) coeff_ar345 = runif(3,-1,1) stazRoot12 = abs(polyroot(c(1,-coeff_ar12))) stazRoot12345 = abs(polyroot(c(1,-coeff_ar12,-coeff_ar345))) if (sum(stazRoot12>1)==2 & sum(stazRoot12345>1)==5) { ds_coeff_ar25[i,] = c(coeff_ar12,coeff_ar345) }}
# teniamo solo i modelli stazionari (in entrambi i casi)ds_coeff_ar25 = ds_coeff_ar25[which(ds_coeff_ar25$theta1!=0 & ds_coeff_ar25$theta2!=0 & ds_coeff_ar25$theta3!=0 & ds_coeff_ar25$theta4!=0 & ds_coeff_ar25$theta5!=0),]ds_coeff_ar25 = ds_coeff_ar25[1:nrep,]
# lista che conterr? tutti i risultati (sar? lunga nrep)allResultsList = list()
# liste che conterranno, per ognuna delle 4 distribuzioni di errore,# i dati generati dagli AR3 per ogni replicazione (saranno lunghe nsim)yAr25Norm = list()yAr25ShGamma = list()yAr25DoubExp = list()yAr25T = list()
# liste che conterranno, per ognuna delle 4 distribuzioni di errore,# i 5 modelli AR stimati sui dati generati precedentemente (saranno lunghe 5)fitArNorm = list()fitArShGamma = list()fitArDoubExp = list()fitArT = list()
# liste che conterranno, per ognuna delle 4 distribuzioni di errore,# le previsioni fatte sulle ultime 25 osservazioniprevArNorm = list()prevArShGamma = list()prevArDoubExp = list()prevArT = list()
cat("\n")beginTime0 = format(Sys.time(), "%Y-%m-%d %H:%M:%S")beginTime0Math = Sys.time()cat("La procedura inizia in data e ora", beginTime0, "\n")cat("\n")
# ciclo che fa tutto per ognuno degli nrep modelli AR3for (r in 1:nrep) { cat("\n") beginTime = format(Sys.time(), "%Y-%m-%d %H:%M:%S") beginTimeMath = Sys.time() cat("Inizia la replicazione", r, "in data e ora", beginTime, "\n") cat("\n") # ciclo che per ogni simulazione simula i dati, stima i modelli e fa le previsioni (nsim volte) cat("Inizio generazione dei dati, stima dei modelli e calcolo delle previsioni\n") for (k in 1:nsim) { # per ognuna delle 4 distribuzioni generiamo 125 dati dall'AR2/5 con coefficienti ds_coeff_ar25[r,] yAr25Norm[[k]] = ts(c(arima.sim(model=list(order=c(2,0,0), ar=as.vector(t(ds_coeff_ar25[r,1:2]))), n=115, sd=sqrt(sigma2)), arima.sim(model=list(order=c(5,0,0), ar=as.vector(t(ds_coeff_ar25[r,]))), n=10, sd=sqrt(sigma2)))) yAr25ShGamma[[k]] = ts(c(arima.sim(model=list(order=c(2,0,0), ar=as.vector(t(ds_coeff_ar25[r,1:2]))), n=115, rand.gen=function(n,...) sqrt(sigma2/3)*(rgamma(n, shape=3, scale=1) - 3)), arima.sim(model=list(order=c(5,0,0), ar=as.vector(t(ds_coeff_ar25[r,]))), n=10, rand.gen=function(n,...) sqrt(sigma2/3)*(rgamma(n, shape=3, scale=1) - 3)))) yAr25DoubExp[[k]] = ts(c(arima.sim(model=list(order=c(2,0,0), ar=as.vector(t(ds_coeff_ar25[r,1:2]))), n=115, rand.gen=function(n,...) sqrt(sigma2/2)*rlaplace(n, location=0, scale=1)), arima.sim(model=list(order=c(5,0,0), ar=as.vector(t(ds_coeff_ar25[r,]))), n=10, rand.gen=function(n,...) sqrt(sigma2/2)*rlaplace(n, location=0, scale=1)))) yAr25T[[k]] = ts(c(arima.sim(model=list(order=c(2,0,0), ar=as.vector(t(ds_coeff_ar25[r,1:2]))), n=115, rand.gen=function(n,...) sqrt(sigma2/2)*rt(n, df=4)), arima.sim(model=list(order=c(5,0,0), ar=as.vector(t(ds_coeff_ar25[r,]))), n=10, rand.gen=function(n,...) sqrt(sigma2/2)*rt(n, df=4)))) # stimiamo i 5 modelli AR sui primi 100 dati e con questi facciamo le previsioni sugli ultimi 25 dati # generati dall'AR2/5 con errori normali # AR1Norm # proviamo il metodo CSS-ML AR1Norm = try(arima(x=yAr25Norm[[k]][1:100],order=c(1,0,0), method="CSS-ML"), silent=FALSE) # se CSS-ML d? errore, proviamo il metodo ML if (class(AR1Norm)=="try-error") { AR1Norm = try(arima(x=yAr25Norm[[k]][1:100],order=c(1,0,0), method="ML"), silent=FALSE) # se anche ML d? errore, proviamo il metodo CSS if (class(AR1Norm)=="try-error") { AR1Norm = try(arima(x=yAr25Norm[[k]][1:100],order=c(1,0,0), method="CSS"), silent=FALSE) } # se ML non d? errore, controlliamo le varianze: se ce n'?? almeno una negativa, proviamo il metodo CSS else { if (sum(diag(AR1Norm$var.coef)>0)<2) { AR1Norm = try(arima(x=yAr25Norm[[k]][1:100],order=c(1,0,0), method="CSS"), silent=FALSE) } } } # se CSS-ML non d? errore, controlliamo le varianze: se ce n'?? almeno una negativa, proviamo il metodo ML else { if (sum(diag(AR1Norm$var.coef)>0)<2) { AR1Norm = try(arima(x=yAr25Norm[[k]][1:100],order=c(1,0,0), method="ML"), silent=FALSE) # se anche ML d? errore, proviamo il metodo CSS if (class(AR1Norm)=="try-error") { AR1Norm = try(arima(x=yAr25Norm[[k]][1:100],order=c(1,0,0), method="CSS"), silent=FALSE) } # se ML non d? errore, controlliamo le varianze: se ce n'?? almeno una negativa, proviamo il metodo CSS else { if (sum(diag(AR1Norm$var.coef)>0)<2) { AR1Norm = try(arima(x=yAr25Norm[[k]][1:100],order=c(1,0,0), method="CSS"), silent=FALSE) } } } } # AR2Norm # proviamo il metodo CSS-ML AR2Norm = try(arima(x=yAr25Norm[[k]][1:100],order=c(2,0,0), method="CSS-ML"), silent=FALSE) # se CSS-ML d? errore, proviamo il metodo ML if (class(AR2Norm)=="try-error") { AR2Norm = try(arima(x=yAr25Norm[[k]][1:100],order=c(2,0,0), method="ML"), silent=FALSE) # se anche ML d? errore, proviamo il metodo CSS if (class(AR2Norm)=="try-error") { AR2Norm = try(arima(x=yAr25Norm[[k]][1:100],order=c(2,0,0), method="CSS"), silent=FALSE) } # se ML non d? errore, controlliamo le varianze: se ce n'?? almeno una negativa, proviamo il metodo CSS else { if (sum(diag(AR2Norm$var.coef)>0)<2) { AR2Norm = try(arima(x=yAr25Norm[[k]][1:100],order=c(2,0,0), method="CSS"), silent=FALSE) } } } # se CSS-ML non d? errore, controlliamo le varianze: se ce n'?? almeno una negativa, proviamo il metodo ML else { if (sum(diag(AR2Norm$var.coef)>0)<2) { AR2Norm = try(arima(x=yAr25Norm[[k]][1:100],order=c(2,0,0), method="ML"), silent=FALSE) # se anche ML d? errore, proviamo il metodo CSS if (class(AR2Norm)=="try-error") { AR2Norm = try(arima(x=yAr25Norm[[k]][1:100],order=c(2,0,0), method="CSS"), silent=FALSE) } # se ML non d? errore, controlliamo le varianze: se ce n'?? almeno una negativa, proviamo il metodo CSS else { if (sum(diag(AR2Norm$var.coef)>0)<2) { AR2Norm = try(arima(x=yAr25Norm[[k]][1:100],order=c(2,0,0), method="CSS"), silent=FALSE) } } } } # AR3Norm # proviamo il metodo CSS-ML AR3Norm = try(arima(x=yAr25Norm[[k]][1:100],order=c(3,0,0), method="CSS-ML"), silent=FALSE) # se CSS-ML d? errore, proviamo il metodo ML if (class(AR3Norm)=="try-error") { AR3Norm = try(arima(x=yAr25Norm[[k]][1:100],order=c(3,0,0), method="ML"), silent=FALSE) # se anche ML d? errore, proviamo il metodo CSS if (class(AR3Norm)=="try-error") { AR3Norm = try(arima(x=yAr25Norm[[k]][1:100],order=c(3,0,0), method="CSS"), silent=FALSE) } # se ML non d? errore, controlliamo le varianze: se ce n'?? almeno una negativa, proviamo il metodo CSS else { if (sum(diag(AR3Norm$var.coef)>0)<2) { AR3Norm = try(arima(x=yAr25Norm[[k]][1:100],order=c(3,0,0), method="CSS"), silent=FALSE) } } } # se CSS-ML non d? errore, controlliamo le varianze: se ce n'?? almeno una negativa, proviamo il metodo ML else { if (sum(diag(AR3Norm$var.coef)>0)<2) { AR3Norm = try(arima(x=yAr25Norm[[k]][1:100],order=c(3,0,0), method="ML"), silent=FALSE) # se anche ML d? errore, proviamo il metodo CSS if (class(AR3Norm)=="try-error") { AR3Norm = try(arima(x=yAr25Norm[[k]][1:100],order=c(3,0,0), method="CSS"), silent=FALSE) } # se ML non d? errore, controlliamo le varianze: se ce n'?? almeno una negativa, proviamo il metodo CSS else { if (sum(diag(AR3Norm$var.coef)>0)<2) { AR3Norm = try(arima(x=yAr25Norm[[k]][1:100],order=c(3,0,0), method="CSS"), silent=FALSE) } } } } # AR4Norm # proviamo il metodo CSS-ML AR4Norm = try(arima(x=yAr25Norm[[k]][1:100],order=c(4,0,0), method="CSS-ML"), silent=FALSE) # se CSS-ML d? errore, proviamo il metodo ML if (class(AR4Norm)=="try-error") { AR4Norm = try(arima(x=yAr25Norm[[k]][1:100],order=c(4,0,0), method="ML"), silent=FALSE) # se anche ML d? errore, proviamo il metodo CSS if (class(AR4Norm)=="try-error") { AR4Norm = try(arima(x=yAr25Norm[[k]][1:100],order=c(4,0,0), method="CSS"), silent=FALSE) } # se ML non d? errore, controlliamo le varianze: se ce n'?? almeno una negativa, proviamo il metodo CSS else { if (sum(diag(AR4Norm$var.coef)>0)<2) { AR4Norm = try(arima(x=yAr25Norm[[k]][1:100],order=c(4,0,0), method="CSS"), silent=FALSE) } } } # se CSS-ML non d? errore, controlliamo le varianze: se ce n'?? almeno una negativa, proviamo il metodo ML else { if (sum(diag(AR4Norm$var.coef)>0)<2) { AR4Norm = try(arima(x=yAr25Norm[[k]][1:100],order=c(4,0,0), method="ML"), silent=FALSE) # se anche ML d? errore, proviamo il metodo CSS if (class(AR4Norm)=="try-error") { AR4Norm = try(arima(x=yAr25Norm[[k]][1:100],order=c(4,0,0), method="CSS"), silent=FALSE) } # se ML non d? errore, controlliamo le varianze: se ce n'?? almeno una negativa, proviamo il metodo CSS else { if (sum(diag(AR4Norm$var.coef)>0)<2) { AR4Norm = try(arima(x=yAr25Norm[[k]][1:100],order=c(4,0,0), method="CSS"), silent=FALSE) } } } } # AR5Norm # proviamo il metodo CSS-ML AR5Norm = try(arima(x=yAr25Norm[[k]][1:100],order=c(5,0,0), method="CSS-ML"), silent=FALSE) # se CSS-ML d? errore, proviamo il metodo ML if (class(AR5Norm)=="try-error") { AR5Norm = try(arima(x=yAr25Norm[[k]][1:100],order=c(5,0,0), method="ML"), silent=FALSE) # se anche ML d? errore, proviamo il metodo CSS if (class(AR5Norm)=="try-error") { AR5Norm = try(arima(x=yAr25Norm[[k]][1:100],order=c(5,0,0), method="CSS"), silent=FALSE) } # se ML non d? errore, controlliamo le varianze: se ce n'?? almeno una negativa, proviamo il metodo CSS else { if (sum(diag(AR5Norm$var.coef)>0)<2) { AR5Norm = try(arima(x=yAr25Norm[[k]][1:100],order=c(5,0,0), method="CSS"), silent=FALSE) } } } # se CSS-ML non d? errore, controlliamo le varianze: se ce n'?? almeno una negativa, proviamo il metodo ML else { if (sum(diag(AR5Norm$var.coef)>0)<2) { AR5Norm = try(arima(x=yAr25Norm[[k]][1:100],order=c(5,0,0), method="ML"), silent=FALSE) # se anche ML d? errore, proviamo il metodo CSS if (class(AR5Norm)=="try-error") { AR5Norm = try(arima(x=yAr25Norm[[k]][1:100],order=c(5,0,0), method="CSS"), silent=FALSE) } # se ML non d? errore, controlliamo le varianze: se ce n'?? almeno una negativa, proviamo il metodo CSS else { if (sum(diag(AR5Norm$var.coef)>0)<2) { AR5Norm = try(arima(x=yAr25Norm[[k]][1:100],order=c(5,0,0), method="CSS"), silent=FALSE) } } } } fitArNorm[[k]] = list(AR1 = AR1Norm, AR2 = AR2Norm, AR3 = AR3Norm, AR4 = AR4Norm, AR5 = AR5Norm) prevArNorm[[k]] = list(AR1 = predict(object=fitArNorm[[k]]$AR1, n.ahead=25), AR2 = predict(object=fitArNorm[[k]]$AR2, n.ahead=25), AR3 = predict(object=fitArNorm[[k]]$AR3, n.ahead=25), AR4 = predict(object=fitArNorm[[k]]$AR4, n.ahead=25), AR5 = predict(object=fitArNorm[[k]]$AR5, n.ahead=25) ) # stimiamo i 5 modelli AR sui primi 100 dati e con questi facciamo le previsioni sugli ultimi 25 dati # generati dall'AR3 con errori gamma traslata # AR1ShGamma # proviamo il metodo CSS-ML AR1ShGamma = try(arima(x=yAr25ShGamma[[k]][1:100],order=c(1,0,0), method="CSS-ML"), silent=FALSE) # se CSS-ML d? errore, proviamo il metodo ML if (class(AR1ShGamma)=="try-error") { AR1ShGamma = try(arima(x=yAr25ShGamma[[k]][1:100],order=c(1,0,0), method="ML"), silent=FALSE) # se anche ML d? errore, proviamo il metodo CSS if (class(AR1ShGamma)=="try-error") { AR1ShGamma = try(arima(x=yAr25ShGamma[[k]][1:100],order=c(1,0,0), method="CSS"), silent=FALSE) } # se ML non d? errore, controlliamo le varianze: se ce n'?? almeno una negativa, proviamo il metodo CSS else { if (sum(diag(AR1ShGamma$var.coef)>0)<2) { AR1ShGamma = try(arima(x=yAr25ShGamma[[k]][1:100],order=c(1,0,0), method="CSS"), silent=FALSE) } } } # se CSS-ML non d? errore, controlliamo le varianze: se ce n'?? almeno una negativa, proviamo il metodo ML else { if (sum(diag(AR1ShGamma$var.coef)>0)<2) { AR1ShGamma = try(arima(x=yAr25ShGamma[[k]][1:100],order=c(1,0,0), method="ML"), silent=FALSE) # se anche ML d? errore, proviamo il metodo CSS if (class(AR1ShGamma)=="try-error") { AR1ShGamma = try(arima(x=yAr25ShGamma[[k]][1:100],order=c(1,0,0), method="CSS"), silent=FALSE) } # se ML non d? errore, controlliamo le varianze: se ce n'?? almeno una negativa, proviamo il metodo CSS else { if (sum(diag(AR1ShGamma$var.coef)>0)<2) { AR1ShGamma = try(arima(x=yAr25ShGamma[[k]][1:100],order=c(1,0,0), method="CSS"), silent=FALSE) } } } } # AR2ShGamma # proviamo il metodo CSS-ML AR2ShGamma = try(arima(x=yAr25ShGamma[[k]][1:100],order=c(2,0,0), method="CSS-ML"), silent=FALSE) # se CSS-ML d? errore, proviamo il metodo ML if (class(AR2ShGamma)=="try-error") { AR2ShGamma = try(arima(x=yAr25ShGamma[[k]][1:100],order=c(2,0,0), method="ML"), silent=FALSE) # se anche ML d? errore, proviamo il metodo CSS if (class(AR2ShGamma)=="try-error") { AR2ShGamma = try(arima(x=yAr25ShGamma[[k]][1:100],order=c(2,0,0), method="CSS"), silent=FALSE) } # se ML non d? errore, controlliamo le varianze: se ce n'?? almeno una negativa, proviamo il metodo CSS else { if (sum(diag(AR2ShGamma$var.coef)>0)<2) { AR2ShGamma = try(arima(x=yAr25ShGamma[[k]][1:100],order=c(2,0,0), method="CSS"), silent=FALSE) } } } # se CSS-ML non d? errore, controlliamo le varianze: se ce n'?? almeno una negativa, proviamo il metodo ML else { if (sum(diag(AR2ShGamma$var.coef)>0)<2) { AR2ShGamma = try(arima(x=yAr25ShGamma[[k]][1:100],order=c(2,0,0), method="ML"), silent=FALSE) # se anche ML d? errore, proviamo il metodo CSS if (class(AR2ShGamma)=="try-error") { AR2ShGamma = try(arima(x=yAr25ShGamma[[k]][1:100],order=c(2,0,0), method="CSS"), silent=FALSE) } # se ML non d? errore, controlliamo le varianze: se ce n'?? almeno una negativa, proviamo il metodo CSS else { if (sum(diag(AR2ShGamma$var.coef)>0)<2) { AR2ShGamma = try(arima(x=yAr25ShGamma[[k]][1:100],order=c(2,0,0), method="CSS"), silent=FALSE) } } } } # AR3ShGamma # proviamo il metodo CSS-ML AR3ShGamma = try(arima(x=yAr25ShGamma[[k]][1:100],order=c(3,0,0), method="CSS-ML"), silent=FALSE) # se CSS-ML d? errore, proviamo il metodo ML if (class(AR3ShGamma)=="try-error") { AR3ShGamma = try(arima(x=yAr25ShGamma[[k]][1:100],order=c(3,0,0), method="ML"), silent=FALSE) # se anche ML d? errore, proviamo il metodo CSS if (class(AR3ShGamma)=="try-error") { AR3ShGamma = try(arima(x=yAr25ShGamma[[k]][1:100],order=c(3,0,0), method="CSS"), silent=FALSE) } # se ML non d? errore, controlliamo le varianze: se ce n'?? almeno una negativa, proviamo il metodo CSS else { if (sum(diag(AR3ShGamma$var.coef)>0)<2) { AR3ShGamma = try(arima(x=yAr25ShGamma[[k]][1:100],order=c(3,0,0), method="CSS"), silent=FALSE) } } } # se CSS-ML non d? errore, controlliamo le varianze: se ce n'?? almeno una negativa, proviamo il metodo ML else { if (sum(diag(AR3ShGamma$var.coef)>0)<2) { AR3ShGamma = try(arima(x=yAr25ShGamma[[k]][1:100],order=c(3,0,0), method="ML"), silent=FALSE) # se anche ML d? errore, proviamo il metodo CSS if (class(AR3ShGamma)=="try-error") { AR3ShGamma = try(arima(x=yAr25ShGamma[[k]][1:100],order=c(3,0,0), method="CSS"), silent=FALSE) } # se ML non d? errore, controlliamo le varianze: se ce n'?? almeno una negativa, proviamo il metodo CSS else { if (sum(diag(AR3ShGamma$var.coef)>0)<2) { AR3ShGamma = try(arima(x=yAr25ShGamma[[k]][1:100],order=c(3,0,0), method="CSS"), silent=FALSE) } } } } # AR4ShGamma # proviamo il metodo CSS-ML AR4ShGamma = try(arima(x=yAr25ShGamma[[k]][1:100],order=c(4,0,0), method="CSS-ML"), silent=FALSE) # se CSS-ML d? errore, proviamo il metodo ML if (class(AR4ShGamma)=="try-error") { AR4ShGamma = try(arima(x=yAr25ShGamma[[k]][1:100],order=c(4,0,0), method="ML"), silent=FALSE) # se anche ML d? errore, proviamo il metodo CSS if (class(AR4ShGamma)=="try-error") { AR4ShGamma = try(arima(x=yAr25ShGamma[[k]][1:100],order=c(4,0,0), method="CSS"), silent=FALSE) } # se ML non d? errore, controlliamo le varianze: se ce n'?? almeno una negativa, proviamo il metodo CSS else { if (sum(diag(AR4ShGamma$var.coef)>0)<2) { AR4ShGamma = try(arima(x=yAr25ShGamma[[k]][1:100],order=c(4,0,0), method="CSS"), silent=FALSE) } } } # se CSS-ML non d? errore, controlliamo le varianze: se ce n'?? almeno una negativa, proviamo il metodo ML else { if (sum(diag(AR4ShGamma$var.coef)>0)<2) { AR4ShGamma = try(arima(x=yAr25ShGamma[[k]][1:100],order=c(4,0,0), method="ML"), silent=FALSE) # se anche ML d? errore, proviamo il metodo CSS if (class(AR4ShGamma)=="try-error") { AR4ShGamma = try(arima(x=yAr25ShGamma[[k]][1:100],order=c(4,0,0), method="CSS"), silent=FALSE) } # se ML non d? errore, controlliamo le varianze: se ce n'?? almeno una negativa, proviamo il metodo CSS else { if (sum(diag(AR4ShGamma$var.coef)>0)<2) { AR4ShGamma = try(arima(x=yAr25ShGamma[[k]][1:100],order=c(4,0,0), method="CSS"), silent=FALSE) } } } } # AR5ShGamma # proviamo il metodo CSS-ML AR5ShGamma = try(arima(x=yAr25ShGamma[[k]][1:100],order=c(5,0,0), method="CSS-ML"), silent=FALSE) # se CSS-ML d? errore, proviamo il metodo ML if (class(AR5ShGamma)=="try-error") { AR5ShGamma = try(arima(x=yAr25ShGamma[[k]][1:100],order=c(5,0,0), method="ML"), silent=FALSE) # se anche ML d? errore, proviamo il metodo CSS if (class(AR5ShGamma)=="try-error") { AR5ShGamma = try(arima(x=yAr25ShGamma[[k]][1:100],order=c(5,0,0), method="CSS"), silent=FALSE) } # se ML non d? errore, controlliamo le varianze: se ce n'?? almeno una negativa, proviamo il metodo CSS else { if (sum(diag(AR5ShGamma$var.coef)>0)<2) { AR5ShGamma = try(arima(x=yAr25ShGamma[[k]][1:100],order=c(5,0,0), method="CSS"), silent=FALSE) } } } # se CSS-ML non d? errore, controlliamo le varianze: se ce n'?? almeno una negativa, proviamo il metodo ML else { if (sum(diag(AR5ShGamma$var.coef)>0)<2) { AR5ShGamma = try(arima(x=yAr25ShGamma[[k]][1:100],order=c(5,0,0), method="ML"), silent=FALSE) # se anche ML d? errore, proviamo il metodo CSS if (class(AR5ShGamma)=="try-error") { AR5ShGamma = try(arima(x=yAr25ShGamma[[k]][1:100],order=c(5,0,0), method="CSS"), silent=FALSE) } # se ML non d? errore, controlliamo le varianze: se ce n'?? almeno una negativa, proviamo il metodo CSS else { if (sum(diag(AR5ShGamma$var.coef)>0)<2) { AR5ShGamma = try(arima(x=yAr25ShGamma[[k]][1:100],order=c(5,0,0), method="CSS"), silent=FALSE) } } } } fitArShGamma[[k]] = list(AR1 = AR1ShGamma, AR2 = AR2ShGamma, AR3 = AR3ShGamma, AR4 = AR4ShGamma, AR5 = AR5ShGamma) prevArShGamma[[k]] = list(AR1 = predict(object=fitArShGamma[[k]]$AR1, n.ahead=25), AR2 = predict(object=fitArShGamma[[k]]$AR2, n.ahead=25), AR3 = predict(object=fitArShGamma[[k]]$AR3, n.ahead=25), AR4 = predict(object=fitArShGamma[[k]]$AR4, n.ahead=25), AR5 = predict(object=fitArShGamma[[k]]$AR5, n.ahead=25) ) # stimiamo i 5 modelli AR sui primi 100 dati e con questi facciamo le previsioni sugli ultimi 25 dati # generati dall'AR3 con errori Laplace # AR1DoubExp # proviamo il metodo CSS-ML AR1DoubExp = try(arima(x=yAr25DoubExp[[k]][1:100],order=c(1,0,0), method="CSS-ML"), silent=FALSE) # se CSS-ML d? errore, proviamo il metodo ML if (class(AR1DoubExp)=="try-error") { AR1DoubExp = try(arima(x=yAr25DoubExp[[k]][1:100],order=c(1,0,0), method="ML"), silent=FALSE) # se anche ML d? errore, proviamo il metodo CSS if (class(AR1DoubExp)=="try-error") { AR1DoubExp = try(arima(x=yAr25DoubExp[[k]][1:100],order=c(1,0,0), method="CSS"), silent=FALSE) } # se ML non d? errore, controlliamo le varianze: se ce n'?? almeno una negativa, proviamo il metodo CSS else { if (sum(diag(AR1DoubExp$var.coef)>0)<2) { AR1DoubExp = try(arima(x=yAr25DoubExp[[k]][1:100],order=c(1,0,0), method="CSS"), silent=FALSE) } } } # se CSS-ML non d? errore, controlliamo le varianze: se ce n'?? almeno una negativa, proviamo il metodo ML else { if (sum(diag(AR1DoubExp$var.coef)>0)<2) { AR1DoubExp = try(arima(x=yAr25DoubExp[[k]][1:100],order=c(1,0,0), method="ML"), silent=FALSE) # se anche ML d? errore, proviamo il metodo CSS if (class(AR1DoubExp)=="try-error") { AR1DoubExp = try(arima(x=yAr25DoubExp[[k]][1:100],order=c(1,0,0), method="CSS"), silent=FALSE) } # se ML non d? errore, controlliamo le varianze: se ce n'?? almeno una negativa, proviamo il metodo CSS else { if (sum(diag(AR1DoubExp$var.coef)>0)<2) { AR1DoubExp = try(arima(x=yAr25DoubExp[[k]][1:100],order=c(1,0,0), method="CSS"), silent=FALSE) } } } } # AR2DoubExp # proviamo il metodo CSS-ML AR2DoubExp = try(arima(x=yAr25DoubExp[[k]][1:100],order=c(2,0,0), method="CSS-ML"), silent=FALSE) # se CSS-ML d? errore, proviamo il metodo ML if (class(AR2DoubExp)=="try-error") { AR2DoubExp = try(arima(x=yAr25DoubExp[[k]][1:100],order=c(2,0,0), method="ML"), silent=FALSE) # se anche ML d? errore, proviamo il metodo CSS if (class(AR2DoubExp)=="try-error") { AR2DoubExp = try(arima(x=yAr25DoubExp[[k]][1:100],order=c(2,0,0), method="CSS"), silent=FALSE) } # se ML non d? errore, controlliamo le varianze: se ce n'?? almeno una negativa, proviamo il metodo CSS else { if (sum(diag(AR2DoubExp$var.coef)>0)<2) { AR2DoubExp = try(arima(x=yAr25DoubExp[[k]][1:100],order=c(2,0,0), method="CSS"), silent=FALSE) } } } # se CSS-ML non d? errore, controlliamo le varianze: se ce n'?? almeno una negativa, proviamo il metodo ML else { if (sum(diag(AR2DoubExp$var.coef)>0)<2) { AR2DoubExp = try(arima(x=yAr25DoubExp[[k]][1:100],order=c(2,0,0), method="ML"), silent=FALSE) # se anche ML d? errore, proviamo il metodo CSS if (class(AR2DoubExp)=="try-error") { AR2DoubExp = try(arima(x=yAr25DoubExp[[k]][1:100],order=c(2,0,0), method="CSS"), silent=FALSE) } # se ML non d? errore, controlliamo le varianze: se ce n'?? almeno una negativa, proviamo il metodo CSS else { if (sum(diag(AR2DoubExp$var.coef)>0)<2) { AR2DoubExp = try(arima(x=yAr25DoubExp[[k]][1:100],order=c(2,0,0), method="CSS"), silent=FALSE) } } } } # AR3DoubExp # proviamo il metodo CSS-ML AR3DoubExp = try(arima(x=yAr25DoubExp[[k]][1:100],order=c(3,0,0), method="CSS-ML"), silent=FALSE) # se CSS-ML d? errore, proviamo il metodo ML if (class(AR3DoubExp)=="try-error") { AR3DoubExp = try(arima(x=yAr25DoubExp[[k]][1:100],order=c(3,0,0), method="ML"), silent=FALSE) # se anche ML d? errore, proviamo il metodo CSS if (class(AR3DoubExp)=="try-error") { AR3DoubExp = try(arima(x=yAr25DoubExp[[k]][1:100],order=c(3,0,0), method="CSS"), silent=FALSE) } # se ML non d? errore, controlliamo le varianze: se ce n'?? almeno una negativa, proviamo il metodo CSS else { if (sum(diag(AR3DoubExp$var.coef)>0)<2) { AR3DoubExp = try(arima(x=yAr25DoubExp[[k]][1:100],order=c(3,0,0), method="CSS"), silent=FALSE) } } } # se CSS-ML non d? errore, controlliamo le varianze: se ce n'?? almeno una negativa, proviamo il metodo ML else { if (sum(diag(AR3DoubExp$var.coef)>0)<2) { AR3DoubExp = try(arima(x=yAr25DoubExp[[k]][1:100],order=c(3,0,0), method="ML"), silent=FALSE) # se anche ML d? errore, proviamo il metodo CSS if (class(AR3DoubExp)=="try-error") { AR3DoubExp = try(arima(x=yAr25DoubExp[[k]][1:100],order=c(3,0,0), method="CSS"), silent=FALSE) } # se ML non d? errore, controlliamo le varianze: se ce n'?? almeno una negativa, proviamo il metodo CSS else { if (sum(diag(AR3DoubExp$var.coef)>0)<2) { AR3DoubExp = try(arima(x=yAr25DoubExp[[k]][1:100],order=c(3,0,0), method="CSS"), silent=FALSE) } } } } # AR4DoubExp # proviamo il metodo CSS-ML AR4DoubExp = try(arima(x=yAr25DoubExp[[k]][1:100],order=c(4,0,0), method="CSS-ML"), silent=FALSE) # se CSS-ML d? errore, proviamo il metodo ML if (class(AR4DoubExp)=="try-error") { AR4DoubExp = try(arima(x=yAr25DoubExp[[k]][1:100],order=c(4,0,0), method="ML"), silent=FALSE) # se anche ML d? errore, proviamo il metodo CSS if (class(AR4DoubExp)=="try-error") { AR4DoubExp = try(arima(x=yAr25DoubExp[[k]][1:100],order=c(4,0,0), method="CSS"), silent=FALSE) } # se ML non d? errore, controlliamo le varianze: se ce n'?? almeno una negativa, proviamo il metodo CSS else { if (sum(diag(AR4DoubExp$var.coef)>0)<2) { AR4DoubExp = try(arima(x=yAr25DoubExp[[k]][1:100],order=c(4,0,0), method="CSS"), silent=FALSE) } } } # se CSS-ML non d? errore, controlliamo le varianze: se ce n'?? almeno una negativa, proviamo il metodo ML else { if (sum(diag(AR4DoubExp$var.coef)>0)<2) { AR4DoubExp = try(arima(x=yAr25DoubExp[[k]][1:100],order=c(4,0,0), method="ML"), silent=FALSE) # se anche ML d? errore, proviamo il metodo CSS if (class(AR4DoubExp)=="try-error") { AR4DoubExp = try(arima(x=yAr25DoubExp[[k]][1:100],order=c(4,0,0), method="CSS"), silent=FALSE) } # se ML non d? errore, controlliamo le varianze: se ce n'?? almeno una negativa, proviamo il metodo CSS else { if (sum(diag(AR4DoubExp$var.coef)>0)<2) { AR4DoubExp = try(arima(x=yAr25DoubExp[[k]][1:100],order=c(4,0,0), method="CSS"), silent=FALSE) } } } } # AR5DoubExp # proviamo il metodo CSS-ML AR5DoubExp = try(arima(x=yAr25DoubExp[[k]][1:100],order=c(5,0,0), method="CSS-ML"), silent=FALSE) # se CSS-ML d? errore, proviamo il metodo ML if (class(AR5DoubExp)=="try-error") { AR5DoubExp = try(arima(x=yAr25DoubExp[[k]][1:100],order=c(5,0,0), method="ML"), silent=FALSE) # se anche ML d? errore, proviamo il metodo CSS if (class(AR5DoubExp)=="try-error") { AR5DoubExp = try(arima(x=yAr25DoubExp[[k]][1:100],order=c(5,0,0), method="CSS"), silent=FALSE) } # se ML non d? errore, controlliamo le varianze: se ce n'?? almeno una negativa, proviamo il metodo CSS else { if (sum(diag(AR5DoubExp$var.coef)>0)<2) { AR5DoubExp = try(arima(x=yAr25DoubExp[[k]][1:100],order=c(5,0,0), method="CSS"), silent=FALSE) } } } # se CSS-ML non d? errore, controlliamo le varianze: se ce n'?? almeno una negativa, proviamo il metodo ML else { if (sum(diag(AR5DoubExp$var.coef)>0)<2) { AR5DoubExp = try(arima(x=yAr25DoubExp[[k]][1:100],order=c(5,0,0), method="ML"), silent=FALSE) # se anche ML d? errore, proviamo il metodo CSS if (class(AR5DoubExp)=="try-error") { AR5DoubExp = try(arima(x=yAr25DoubExp[[k]][1:100],order=c(5,0,0), method="CSS"), silent=FALSE) } # se ML non d? errore, controlliamo le varianze: se ce n'?? almeno una negativa, proviamo il metodo CSS else { if (sum(diag(AR5DoubExp$var.coef)>0)<2) { AR5DoubExp = try(arima(x=yAr25DoubExp[[k]][1:100],order=c(5,0,0), method="CSS"), silent=FALSE) } } } } fitArDoubExp[[k]] = list(AR1 = AR1DoubExp, AR2 = AR2DoubExp, AR3 = AR3DoubExp, AR4 = AR4DoubExp, AR5 = AR5DoubExp) prevArDoubExp[[k]] = list(AR1 = predict(object=fitArDoubExp[[k]]$AR1, n.ahead=25), AR2 = predict(object=fitArDoubExp[[k]]$AR2, n.ahead=25), AR3 = predict(object=fitArDoubExp[[k]]$AR3, n.ahead=25), AR4 = predict(object=fitArDoubExp[[k]]$AR4, n.ahead=25), AR5 = predict(object=fitArDoubExp[[k]]$AR5, n.ahead=25) ) # stimiamo i 5 modelli AR sui primi 100 dati e con questi facciamo le previsioni sugli ultimi 25 dati # generati dall'AR3 con errori t di Student # AR1T # proviamo il metodo CSS-ML AR1T = try(arima(x=yAr25T[[k]][1:100],order=c(1,0,0), method="CSS-ML"), silent=FALSE) # se CSS-ML d? errore, proviamo il metodo ML if (class(AR1T)=="try-error") { AR1T = try(arima(x=yAr25T[[k]][1:100],order=c(1,0,0), method="ML"), silent=FALSE) # se anche ML d? errore, proviamo il metodo CSS if (class(AR1T)=="try-error") { AR1T = try(arima(x=yAr25T[[k]][1:100],order=c(1,0,0), method="CSS"), silent=FALSE) } # se ML non d? errore, controlliamo le varianze: se ce n'?? almeno una negativa, proviamo il metodo CSS else { if (sum(diag(AR1T$var.coef)>0)<2) { AR1T = try(arima(x=yAr25T[[k]][1:100],order=c(1,0,0), method="CSS"), silent=FALSE) } } } # se CSS-ML non d? errore, controlliamo le varianze: se ce n'?? almeno una negativa, proviamo il metodo ML else { if (sum(diag(AR1T$var.coef)>0)<2) { AR1T = try(arima(x=yAr25T[[k]][1:100],order=c(1,0,0), method="ML"), silent=FALSE) # se anche ML d? errore, proviamo il metodo CSS if (class(AR1T)=="try-error") { AR1T = try(arima(x=yAr25T[[k]][1:100],order=c(1,0,0), method="CSS"), silent=FALSE) } # se ML non d? errore, controlliamo le varianze: se ce n'?? almeno una negativa, proviamo il metodo CSS else { if (sum(diag(AR1T$var.coef)>0)<2) { AR1T = try(arima(x=yAr25T[[k]][1:100],order=c(1,0,0), method="CSS"), silent=FALSE) } } } } # AR2T # proviamo il metodo CSS-ML AR2T = try(arima(x=yAr25T[[k]][1:100],order=c(2,0,0), method="CSS-ML"), silent=FALSE) # se CSS-ML d? errore, proviamo il metodo ML if (class(AR2T)=="try-error") { AR2T = try(arima(x=yAr25T[[k]][1:100],order=c(2,0,0), method="ML"), silent=FALSE) # se anche ML d? errore, proviamo il metodo CSS if (class(AR2T)=="try-error") { AR2T = try(arima(x=yAr25T[[k]][1:100],order=c(2,0,0), method="CSS"), silent=FALSE) } # se ML non d? errore, controlliamo le varianze: se ce n'?? almeno una negativa, proviamo il metodo CSS else { if (sum(diag(AR2T$var.coef)>0)<2) { AR2T = try(arima(x=yAr25T[[k]][1:100],order=c(2,0,0), method="CSS"), silent=FALSE) } } } # se CSS-ML non d? errore, controlliamo le varianze: se ce n'?? almeno una negativa, proviamo il metodo ML else { if (sum(diag(AR2T$var.coef)>0)<2) { AR2T = try(arima(x=yAr25T[[k]][1:100],order=c(2,0,0), method="ML"), silent=FALSE) # se anche ML d? errore, proviamo il metodo CSS if (class(AR2T)=="try-error") { AR2T = try(arima(x=yAr25T[[k]][1:100],order=c(2,0,0), method="CSS"), silent=FALSE) } # se ML non d? errore, controlliamo le varianze: se ce n'?? almeno una negativa, proviamo il metodo CSS else { if (sum(diag(AR2T$var.coef)>0)<2) { AR2T = try(arima(x=yAr25T[[k]][1:100],order=c(2,0,0), method="CSS"), silent=FALSE) } } } } # AR3T # proviamo il metodo CSS-ML AR3T = try(arima(x=yAr25T[[k]][1:100],order=c(3,0,0), method="CSS-ML"), silent=FALSE) # se CSS-ML d? errore, proviamo il metodo ML if (class(AR3T)=="try-error") { AR3T = try(arima(x=yAr25T[[k]][1:100],order=c(3,0,0), method="ML"), silent=FALSE) # se anche ML d? errore, proviamo il metodo CSS if (class(AR3T)=="try-error") { AR3T = try(arima(x=yAr25T[[k]][1:100],order=c(3,0,0), method="CSS"), silent=FALSE) } # se ML non d? errore, controlliamo le varianze: se ce n'?? almeno una negativa, proviamo il metodo CSS else { if (sum(diag(AR3T$var.coef)>0)<2) { AR3T = try(arima(x=yAr25T[[k]][1:100],order=c(3,0,0), method="CSS"), silent=FALSE) } } } # se CSS-ML non d? errore, controlliamo le varianze: se ce n'?? almeno una negativa, proviamo il metodo ML else { if (sum(diag(AR3T$var.coef)>0)<2) { AR3T = try(arima(x=yAr25T[[k]][1:100],order=c(3,0,0), method="ML"), silent=FALSE) # se anche ML d? errore, proviamo il metodo CSS if (class(AR3T)=="try-error") { AR3T = try(arima(x=yAr25T[[k]][1:100],order=c(3,0,0), method="CSS"), silent=FALSE) } # se ML non d? errore, controlliamo le varianze: se ce n'?? almeno una negativa, proviamo il metodo CSS else { if (sum(diag(AR3T$var.coef)>0)<2) { AR3T = try(arima(x=yAr25T[[k]][1:100],order=c(3,0,0), method="CSS"), silent=FALSE) } } } } # AR4T # proviamo il metodo CSS-ML AR4T = try(arima(x=yAr25T[[k]][1:100],order=c(4,0,0), method="CSS-ML"), silent=FALSE) # se CSS-ML d? errore, proviamo il metodo ML if (class(AR4T)=="try-error") { AR4T = try(arima(x=yAr25T[[k]][1:100],order=c(4,0,0), method="ML"), silent=FALSE) # se anche ML d? errore, proviamo il metodo CSS if (class(AR4T)=="try-error") { AR4T = try(arima(x=yAr25T[[k]][1:100],order=c(4,0,0), method="CSS"), silent=FALSE) } # se ML non d? errore, controlliamo le varianze: se ce n'?? almeno una negativa, proviamo il metodo CSS else { if (sum(diag(AR4T$var.coef)>0)<2) { AR4T = try(arima(x=yAr25T[[k]][1:100],order=c(4,0,0), method="CSS"), silent=FALSE) } } } # se CSS-ML non d? errore, controlliamo le varianze: se ce n'?? almeno una negativa, proviamo il metodo ML else { if (sum(diag(AR4T$var.coef)>0)<2) { AR4T = try(arima(x=yAr25T[[k]][1:100],order=c(4,0,0), method="ML"), silent=FALSE) # se anche ML d? errore, proviamo il metodo CSS if (class(AR4T)=="try-error") { AR4T = try(arima(x=yAr25T[[k]][1:100],order=c(4,0,0), method="CSS"), silent=FALSE) } # se ML non d? errore, controlliamo le varianze: se ce n'?? almeno una negativa, proviamo il metodo CSS else { if (sum(diag(AR4T$var.coef)>0)<2) { AR4T = try(arima(x=yAr25T[[k]][1:100],order=c(4,0,0), method="CSS"), silent=FALSE) } } } } # AR5T # proviamo il metodo CSS-ML AR5T = try(arima(x=yAr25T[[k]][1:100],order=c(5,0,0), method="CSS-ML"), silent=FALSE) # se CSS-ML d? errore, proviamo il metodo ML if (class(AR5T)=="try-error") { AR5T = try(arima(x=yAr25T[[k]][1:100],order=c(5,0,0), method="ML"), silent=FALSE) # se anche ML d? errore, proviamo il metodo CSS if (class(AR5T)=="try-error") { AR5T = try(arima(x=yAr25T[[k]][1:100],order=c(5,0,0), method="CSS"), silent=FALSE) } # se ML non d? errore, controlliamo le varianze: se ce n'?? almeno una negativa, proviamo il metodo CSS else { if (sum(diag(AR5T$var.coef)>0)<2) { AR5T = try(arima(x=yAr25T[[k]][1:100],order=c(5,0,0), method="CSS"), silent=FALSE) } } } # se CSS-ML non d? errore, controlliamo le varianze: se ce n'?? almeno una negativa, proviamo il metodo ML else { if (sum(diag(AR5T$var.coef)>0)<2) { AR5T = try(arima(x=yAr25T[[k]][1:100],order=c(5,0,0), method="ML"), silent=FALSE) # se anche ML d? errore, proviamo il metodo CSS if (class(AR5T)=="try-error") { AR5T = try(arima(x=yAr25T[[k]][1:100],order=c(5,0,0), method="CSS"), silent=FALSE) } # se ML non d? errore, controlliamo le varianze: se ce n'?? almeno una negativa, proviamo il metodo CSS else { if (sum(diag(AR5T$var.coef)>0)<2) { AR5T = try(arima(x=yAr25T[[k]][1:100],order=c(5,0,0), method="CSS"), silent=FALSE) } } } } fitArT[[k]] = list(AR1 = AR1T, AR2 = AR2T, AR3 = AR3T, AR4 = AR4T, AR5 = AR5T) prevArT[[k]] = list(AR1 = predict(object=fitArT[[k]]$AR1, n.ahead=25), AR2 = predict(object=fitArT[[k]]$AR2, n.ahead=25), AR3 = predict(object=fitArT[[k]]$AR3, n.ahead=25), AR4 = predict(object=fitArT[[k]]$AR4, n.ahead=25), AR5 = predict(object=fitArT[[k]]$AR5, n.ahead=25) ) if (k %% 10 == 0) print(k) } cat("Fine generazione dei dati, stima dei modelli e calcolo delle previsioni\n") # liste che, per ognuna delle 4 distribuzioni, conterranno i pesi delle previsioni combinate weightListNorm = list() weightListShGamma = list() weightListDoubExp = list() weightListT = list() # numero di previsioni e numero di modelli AR usati per la previsione combinata weightT = 25 M = 5 # ciclo che per ogni simulazione calcoler? i pesi delle previsioni combinate con i 3 metodi cat("Inizio calcolo dei pesi delle previsioni combinate con i 3 metodi\n") for (k in 1:nsim) { # lista che conterr? i pesi nel caso normale weightListNorm[[k]] = list( AFTER = data.frame(AR1=double(weightT), AR2=double(weightT), AR3=double(weightT), AR4=double(weightT), AR5=double(weightT)), L1.AFTER = data.frame(AR1=double(weightT), AR2=double(weightT), AR3=double(weightT), AR4=double(weightT), AR5=double(weightT)), H.AFTER = data.frame(AR1=double(weightT), AR2=double(weightT), AR3=double(weightT), AR4=double(weightT), AR5=double(weightT)) ) weightListNorm[[k]]$AFTER[1,] = rep(1/M,M) weightListNorm[[k]]$L1.AFTER[1,] = rep(1/M,M) weightListNorm[[k]]$H.AFTER[1,] = rep(1/M,M) for (i in 2:weightT) { for (j in 1:M) { weightListNorm[[k]]$L1.AFTER[i,j] = weightListNorm[[k]]$L1.AFTER[i-1,1]/sd(yAr25Norm[[k]][1:(100+i-1)])*exp(-lambda* abs(yAr25Norm[[k]][100+i-1]-prevArNorm[[k]][[j]]$pred[i-1]) /( sd(yAr25Norm[[k]][1:(100+i-1)]))) weightListNorm[[k]]$ AFTER[i,j] = weightListNorm[[k]]$AFTER[i-1,1] /sd(yAr25Norm[[k]][1:(100+i-1)])*exp(- (yAr25Norm[[k]][100+i-1]-prevArNorm[[k]][[j]]$pred[i-1])^2/( 2*var(yAr25Norm[[k]][1:(100+i-1)]))) weightListNorm[[k]]$H.AFTER[i,j] = weightListNorm[[k]]$H.AFTER[i-1,1] /sd(yAr25Norm[[k]][1:(100+i-1)])*exp(-lambda*phis(x=(yAr25Norm[[k]][100+i-1]-prevArNorm[[k]][[j]]$pred[i-1]) /(sqrt( 2)*sd(yAr25Norm[[k]][1:(100+i-1)])), s=s)) } } sumWeightsNorm = list(AFTER = apply(weightListNorm[[k]]$AFTER, 1, sum), L1.AFTER = apply(weightListNorm[[k]]$L1.AFTER, 1, sum), H.AFTER = apply(weightListNorm[[k]]$H.AFTER, 1, sum) ) for (i in 1:weightT) { for (j in 1:M) { weightListNorm[[k]]$AFTER[i,j] = weightListNorm[[k]]$AFTER[i,j] /sumWeightsNorm$AFTER[i] weightListNorm[[k]]$L1.AFTER[i,j] = weightListNorm[[k]]$L1.AFTER[i,j]/sumWeightsNorm$L1.AFTER[i] weightListNorm[[k]]$H.AFTER[i,j] = weightListNorm[[k]]$H.AFTER[i,j] /sumWeightsNorm$H.AFTER[i] } } # lista che conterr? i pesi nel caso gamma traslata weightListShGamma[[k]] = list( AFTER = data.frame(AR1=double(weightT), AR2=double(weightT), AR3=double(weightT), AR4=double(weightT), AR5=double(weightT)), L1.AFTER = data.frame(AR1=double(weightT), AR2=double(weightT), AR3=double(weightT), AR4=double(weightT), AR5=double(weightT)), H.AFTER = data.frame(AR1=double(weightT), AR2=double(weightT), AR3=double(weightT), AR4=double(weightT), AR5=double(weightT)) ) weightListShGamma[[k]]$AFTER[1,] = rep(1/M,M) weightListShGamma[[k]]$L1.AFTER[1,] = rep(1/M,M) weightListShGamma[[k]]$H.AFTER[1,] = rep(1/M,M) for (i in 2:weightT) { for (j in 1:M) { weightListShGamma[[k]]$L1.AFTER[i,j] = weightListShGamma[[k]]$L1.AFTER[i-1,1]/sd(yAr25ShGamma[[k]][1:(100+i-1)])*exp(-lambda* abs(yAr25ShGamma[[k]][100+i-1]-prevArShGamma[[k]][[j]]$pred[i-1]) /( sd(yAr25ShGamma[[k]][1:(100+i-1)]))) weightListShGamma[[k]]$ AFTER[i,j] = weightListShGamma[[k]]$AFTER[i-1,1] /sd(yAr25ShGamma[[k]][1:(100+i-1)])*exp(- (yAr25ShGamma[[k]][100+i-1]-prevArShGamma[[k]][[j]]$pred[i-1])^2/( 2*var(yAr25ShGamma[[k]][1:(100+i-1)]))) weightListShGamma[[k]]$H.AFTER[i,j] = weightListShGamma[[k]]$H.AFTER[i-1,1] /sd(yAr25ShGamma[[k]][1:(100+i-1)])*exp(-lambda*phis(x=(yAr25ShGamma[[k]][100+i-1]-prevArShGamma[[k]][[j]]$pred[i-1]) /(sqrt( 2)*sd(yAr25ShGamma[[k]][1:(100+i-1)])), s=s)) } } sumWeightsShGamma = list(AFTER = apply(weightListShGamma[[k]]$AFTER, 1, sum), L1.AFTER = apply(weightListShGamma[[k]]$L1.AFTER, 1, sum), H.AFTER = apply(weightListShGamma[[k]]$H.AFTER, 1, sum) ) for (i in 1:weightT) { for (j in 1:M) { weightListShGamma[[k]]$AFTER[i,j] = weightListShGamma[[k]]$AFTER[i,j] /sumWeightsShGamma$AFTER[i] weightListShGamma[[k]]$L1.AFTER[i,j] = weightListShGamma[[k]]$L1.AFTER[i,j]/sumWeightsShGamma$L1.AFTER[i] weightListShGamma[[k]]$H.AFTER[i,j] = weightListShGamma[[k]]$H.AFTER[i,j] /sumWeightsShGamma$H.AFTER[i] } } # lista che conterr? i pesi nel caso Laplace weightListDoubExp[[k]] = list( AFTER = data.frame(AR1=double(weightT), AR2=double(weightT), AR3=double(weightT), AR4=double(weightT), AR5=double(weightT)), L1.AFTER = data.frame(AR1=double(weightT), AR2=double(weightT), AR3=double(weightT), AR4=double(weightT), AR5=double(weightT)), H.AFTER = data.frame(AR1=double(weightT), AR2=double(weightT), AR3=double(weightT), AR4=double(weightT), AR5=double(weightT)) ) weightListDoubExp[[k]]$AFTER[1,] = rep(1/M,M) weightListDoubExp[[k]]$L1.AFTER[1,] = rep(1/M,M) weightListDoubExp[[k]]$H.AFTER[1,] = rep(1/M,M) for (i in 2:weightT) { for (j in 1:M) { weightListDoubExp[[k]]$L1.AFTER[i,j] = weightListDoubExp[[k]]$L1.AFTER[i-1,1]/sd(yAr25DoubExp[[k]][1:(100+i-1)])*exp(-lambda* abs(yAr25DoubExp[[k]][100+i-1]-prevArDoubExp[[k]][[j]]$pred[i-1]) /( sd(yAr25DoubExp[[k]][1:(100+i-1)]))) weightListDoubExp[[k]]$ AFTER[i,j] = weightListDoubExp[[k]]$AFTER[i-1,1] /sd(yAr25DoubExp[[k]][1:(100+i-1)])*exp(- (yAr25DoubExp[[k]][100+i-1]-prevArDoubExp[[k]][[j]]$pred[i-1])^2/( 2*var(yAr25DoubExp[[k]][1:(100+i-1)]))) weightListDoubExp[[k]]$H.AFTER[i,j] = weightListDoubExp[[k]]$H.AFTER[i-1,1] /sd(yAr25DoubExp[[k]][1:(100+i-1)])*exp(-lambda*phis(x=(yAr25DoubExp[[k]][100+i-1]-prevArDoubExp[[k]][[j]]$pred[i-1]) /(sqrt( 2)*sd(yAr25DoubExp[[k]][1:(100+i-1)])), s=s)) } } sumWeightsDoubExp = list(AFTER = apply(weightListDoubExp[[k]]$AFTER, 1, sum), L1.AFTER = apply(weightListDoubExp[[k]]$L1.AFTER, 1, sum), H.AFTER = apply(weightListDoubExp[[k]]$H.AFTER, 1, sum) ) for (i in 1:weightT) { for (j in 1:M) { weightListDoubExp[[k]]$AFTER[i,j] = weightListDoubExp[[k]]$AFTER[i,j] /sumWeightsDoubExp$AFTER[i] weightListDoubExp[[k]]$L1.AFTER[i,j] = weightListDoubExp[[k]]$L1.AFTER[i,j]/sumWeightsDoubExp$L1.AFTER[i] weightListDoubExp[[k]]$H.AFTER[i,j] = weightListDoubExp[[k]]$H.AFTER[i,j] /sumWeightsDoubExp$H.AFTER[i] } } # lista che conterr? i pesi nel caso t di Student weightListT[[k]] = list( AFTER = data.frame(AR1=double(weightT), AR2=double(weightT), AR3=double(weightT), AR4=double(weightT), AR5=double(weightT)), L1.AFTER = data.frame(AR1=double(weightT), AR2=double(weightT), AR3=double(weightT), AR4=double(weightT), AR5=double(weightT)), H.AFTER = data.frame(AR1=double(weightT), AR2=double(weightT), AR3=double(weightT), AR4=double(weightT), AR5=double(weightT)) ) weightListT[[k]]$AFTER[1,] = rep(1/M,M) weightListT[[k]]$L1.AFTER[1,] = rep(1/M,M) weightListT[[k]]$H.AFTER[1,] = rep(1/M,M) for (i in 2:weightT) { for (j in 1:M) { weightListT[[k]]$L1.AFTER[i,j] = weightListT[[k]]$L1.AFTER[i-1,1]/sd(yAr25T[[k]][1:(100+i-1)])*exp(-lambda* abs(yAr25T[[k]][100+i-1]-prevArT[[k]][[j]]$pred[i-1]) /( sd(yAr25T[[k]][1:(100+i-1)]))) weightListT[[k]]$ AFTER[i,j] = weightListT[[k]]$AFTER[i-1,1] /sd(yAr25T[[k]][1:(100+i-1)])*exp(- (yAr25T[[k]][100+i-1]-prevArT[[k]][[j]]$pred[i-1])^2/( 2*var(yAr25T[[k]][1:(100+i-1)]))) weightListT[[k]]$H.AFTER[i,j] = weightListT[[k]]$H.AFTER[i-1,1] /sd(yAr25T[[k]][1:(100+i-1)])*exp(-lambda*phis(x=(yAr25T[[k]][100+i-1]-prevArT[[k]][[j]]$pred[i-1]) /(sqrt( 2)*sd(yAr25T[[k]][1:(100+i-1)])), s=s)) } } sumWeightsT = list(AFTER = apply(weightListT[[k]]$AFTER, 1, sum), L1.AFTER = apply(weightListT[[k]]$L1.AFTER, 1, sum), H.AFTER = apply(weightListT[[k]]$H.AFTER, 1, sum) ) for (i in 1:weightT) { for (j in 1:M) { weightListT[[k]]$AFTER[i,j] = weightListT[[k]]$AFTER[i,j] /sumWeightsT$AFTER[i] weightListT[[k]]$L1.AFTER[i,j] = weightListT[[k]]$L1.AFTER[i,j]/sumWeightsT$L1.AFTER[i] weightListT[[k]]$H.AFTER[i,j] = weightListT[[k]]$H.AFTER[i,j] /sumWeightsT$H.AFTER[i] } } if (k %% 10 == 0) print(k) } cat("Fine calcolo dei pesi delle previsioni combinate con i 3 metodi\n") # liste che, per ognuna delle 4 distribuzioni, conterranno le previsioni combinate combinedPredictNorm = list() combinedPredictShGamma = list() combinedPredictDoubExp = list() combinedPredictT = list() cat("Inizio calcolo delle previsioni combinate con i 3 metodi\n") for (k in 1:nsim) { # previsioni combinate per il caso normale combinedPredictNorm[[k]] = list(AFTER = double(weightT), L1.AFTER = double(weightT), H.AFTER = double(weightT)) for (i in 1:weightT) { prevArNormki = c(prevArNorm[[k]]$AR1$pred[i], prevArNorm[[k]]$AR2$pred[i], prevArNorm[[k]]$AR3$pred[i], prevArNorm[[k]]$AR4$pred[i], prevArNorm[[k]]$AR5$pred[i]) combinedPredictNorm[[k]]$AFTER[i] = sum(weightListNorm[[k]]$AFTER[i,] *prevArNormki) combinedPredictNorm[[k]]$L1.AFTER[i] = sum(weightListNorm[[k]]$L1.AFTER[i,]*prevArNormki) combinedPredictNorm[[k]]$H.AFTER[i] = sum(weightListNorm[[k]]$H.AFTER[i,] *prevArNormki) } # previsioni combinate per il caso gamma traslata combinedPredictShGamma[[k]] = list(AFTER = double(weightT), L1.AFTER = double(weightT), H.AFTER = double(weightT)) for (i in 1:weightT) { prevArShGammaki = c(prevArShGamma[[k]]$AR1$pred[i], prevArShGamma[[k]]$AR2$pred[i], prevArShGamma[[k]]$AR3$pred[i], prevArShGamma[[k]]$AR4$pred[i], prevArShGamma[[k]]$AR5$pred[i]) combinedPredictShGamma[[k]]$AFTER[i] = sum(weightListShGamma[[k]]$AFTER[i,] *prevArShGammaki) combinedPredictShGamma[[k]]$L1.AFTER[i] = sum(weightListShGamma[[k]]$L1.AFTER[i,]*prevArShGammaki) combinedPredictShGamma[[k]]$H.AFTER[i] = sum(weightListShGamma[[k]]$H.AFTER[i,] *prevArShGammaki) } # previsioni combinate per il caso Laplace combinedPredictDoubExp[[k]] = list(AFTER = double(weightT), L1.AFTER = double(weightT), H.AFTER = double(weightT)) for (i in 1:weightT) { prevArDoubExpki = c(prevArDoubExp[[k]]$AR1$pred[i], prevArDoubExp[[k]]$AR2$pred[i], prevArDoubExp[[k]]$AR3$pred[i], prevArDoubExp[[k]]$AR4$pred[i], prevArDoubExp[[k]]$AR5$pred[i]) combinedPredictDoubExp[[k]]$AFTER[i] = sum(weightListDoubExp[[k]]$AFTER[i,] *prevArDoubExpki) combinedPredictDoubExp[[k]]$L1.AFTER[i] = sum(weightListDoubExp[[k]]$L1.AFTER[i,]*prevArDoubExpki) combinedPredictDoubExp[[k]]$H.AFTER[i] = sum(weightListDoubExp[[k]]$H.AFTER[i,] *prevArDoubExpki) } # previsioni combinate per il caso t di Student combinedPredictT[[k]] = list(AFTER = double(weightT), L1.AFTER = double(weightT), H.AFTER = double(weightT)) for (i in 1:weightT) { prevArTki = c(prevArT[[k]]$AR1$pred[i], prevArT[[k]]$AR2$pred[i], prevArT[[k]]$AR3$pred[i], prevArT[[k]]$AR4$pred[i], prevArT[[k]]$AR5$pred[i]) combinedPredictT[[k]]$AFTER[i] = sum(weightListT[[k]]$AFTER[i,] *prevArTki) combinedPredictT[[k]]$L1.AFTER[i] = sum(weightListT[[k]]$L1.AFTER[i,]*prevArTki) combinedPredictT[[k]]$H.AFTER[i] = sum(weightListT[[k]]$H.AFTER[i,] *prevArTki) } if (k %% 10 == 0) print(k) } cat("Fine calcolo delle previsioni combinate con i 3 metodi\n") # MSE nel caso normale (delle 25 osservazioni previste si usano le ultime 20) cat("Inizio calcolo degli MSE con i 3 metodi\n") MSENorm = data.frame(AFTER = double(nsim), L1.AFTER = double(nsim), H.AFTER = double(nsim)) for (k in 1:nsim) { MSENorm[k,"AFTER"] = mean((yAr25Norm[[k]][106:125]-combinedPredictNorm[[k]]$AFTER[6:25])^2) MSENorm[k,"L1.AFTER"] = mean((yAr25Norm[[k]][106:125]-combinedPredictNorm[[k]]$L1.AFTER[6:25])^2) MSENorm[k,"H.AFTER"] = mean((yAr25Norm[[k]][106:125]-combinedPredictNorm[[k]]$H.AFTER[6:25])^2)# print(k) } # MSE nel caso gamma traslata (delle 25 osservazioni previste si usano le ultime 20) MSEShGamma = data.frame(AFTER = double(nsim), L1.AFTER = double(nsim), H.AFTER = double(nsim)) for (k in 1:nsim) { MSEShGamma[k,"AFTER"] = mean((yAr25ShGamma[[k]][106:125]-combinedPredictShGamma[[k]]$AFTER[6:25])^2) MSEShGamma[k,"L1.AFTER"] = mean((yAr25ShGamma[[k]][106:125]-combinedPredictShGamma[[k]]$L1.AFTER[6:25])^2) MSEShGamma[k,"H.AFTER"] = mean((yAr25ShGamma[[k]][106:125]-combinedPredictShGamma[[k]]$H.AFTER[6:25])^2)# print(k) } # MSE nel caso Laplace (delle 25 osservazioni previste si usano le ultime 20) MSEDoubExp = data.frame(AFTER = double(nsim), L1.AFTER = double(nsim), H.AFTER = double(nsim)) for (k in 1:nsim) { MSEDoubExp[k,"AFTER"] = mean((yAr25DoubExp[[k]][106:125]-combinedPredictDoubExp[[k]]$AFTER[6:25])^2) MSEDoubExp[k,"L1.AFTER"] = mean((yAr25DoubExp[[k]][106:125]-combinedPredictDoubExp[[k]]$L1.AFTER[6:25])^2) MSEDoubExp[k,"H.AFTER"] = mean((yAr25DoubExp[[k]][106:125]-combinedPredictDoubExp[[k]]$H.AFTER[6:25])^2)# print(k) } # MSE nel caso t di Student (delle 25 osservazioni previste si usano le ultime 20) MSET = data.frame(AFTER = double(nsim), L1.AFTER = double(nsim), H.AFTER = double(nsim)) for (k in 1:nsim) { MSET[k,"AFTER"] = mean((yAr25T[[k]][106:125]-combinedPredictT[[k]]$AFTER[6:25])^2) MSET[k,"L1.AFTER"] = mean((yAr25T[[k]][106:125]-combinedPredictT[[k]]$L1.AFTER[6:25])^2) MSET[k,"H.AFTER"] = mean((yAr25T[[k]][106:125]-combinedPredictT[[k]]$H.AFTER[6:25])^2) if (k %% 10 == 0) print(k) } cat("Fine calcolo degli MSE con i 3 metodi\n") # LISTA FINALE, LUNGA nrep, CHE CONTIENE TUTTI I RISULTATI PER OGNI REPLICAZIONE!!! allResultsList[[r]] = list(coeff_ar3 = ds_coeff_ar25[r,], sigma2 = sigma2, lambda = lambda, s = s, yAr25Norm = yAr25Norm, yAr25ShGamma = yAr25ShGamma, yAr25DoubExp = yAr25DoubExp, yAr25T = yAr25T, fitArNorm = fitArNorm, fitArShGamma = fitArShGamma, fitArDoubExp = fitArDoubExp, fitArT = fitArT, prevArNorm = prevArNorm, prevArShGamma = prevArShGamma, prevArDoubExp = prevArDoubExp, prevArT = prevArT, weightListNorm = weightListNorm, weightListShGamma = weightListShGamma, weightListDoubExp = weightListDoubExp, weightListT = weightListT, combinedPredictNorm = combinedPredictNorm, combinedPredictShGamma = combinedPredictShGamma, combinedPredictDoubExp = combinedPredictDoubExp, combinedPredictT = combinedPredictT, MSENorm = MSENorm, MSEShGamma = MSEShGamma, MSEDoubExp = MSEDoubExp, MSET = MSET) # salvataggio dei dataframe con gli MSE, identificati per varianza e replicazione (da 1 a nrep) write.table(allResultsList[[r]]$MSENorm, file=paste0(outputpath, "MSENorm-varianza", sigma2, "-modello", r, ".csv"), sep="|", quote=FALSE, row.names=FALSE) write.table(allResultsList[[r]]$MSEShGamma, file=paste0(outputpath, "MSEShGamma-varianza", sigma2, "-modello", r, ".csv"), sep="|", quote=FALSE, row.names=FALSE) write.table(allResultsList[[r]]$MSEDoubExp, file=paste0(outputpath, "MSEDoubExp-varianza", sigma2, "-modello", r, ".csv"), sep="|", quote=FALSE, row.names=FALSE) write.table(allResultsList[[r]]$MSET, file=paste0(outputpath, "MSET-varianza", sigma2, "-modello", r, ".csv"), sep="|", quote=FALSE, row.names=FALSE) cat("Finisce la replicazione", r, "\n") cat("\n") endTime = format(Sys.time(), "%Y-%m-%d %H:%M:%S") endTimeMath = Sys.time() cat("Finisce la replicazione", r, "in data e ora", endTime, "\n") print(endTimeMath-beginTimeMath) cat("\n")
}
endTime0 = format(Sys.time(), "%Y-%m-%d %H:%M:%S")endTime0Math = Sys.time()cat("La procedura ?? iniziata in data e ora", beginTime0, "\n")cat("La procedura ?? finita in data e ora", endTime0, "\n")print(endTime0Math-beginTime0Math)cat("\n")
# salvataggio del dataset con i coefficienti degli AR3 da cui abbiamo simulatowrite.table(ds_coeff_ar25, file=paste0(outputpath, "coefficientiAR25veriBreak-varianza", sigma2, ".csv"), sep="|", quote=FALSE, row.names=FALSE)
# salvataggio della LISTA FINALE CON TUTTI I RISULTATI per ogni varianza utilizzatasave(allResultsList, file=paste0(outputpath, "risultatiBreak-varianza", sigma2, ".Rdata"))
# # summary di tutto insiemesummary(MSENorm[,"L1.AFTER"]/MSENorm[,"AFTER"])summary(MSENorm[,"H.AFTER"]/MSENorm[,"AFTER"])summary(MSEShGamma[,"L1.AFTER"]/MSEShGamma[,"AFTER"])summary(MSEShGamma[,"H.AFTER"]/MSEShGamma[,"AFTER"])summary(MSEDoubExp[,"L1.AFTER"]/MSEDoubExp[,"AFTER"])summary(MSEDoubExp[,"H.AFTER"]/MSEDoubExp[,"AFTER"])summary(MSET[,"L1.AFTER"]/MSET[,"AFTER"])summary(MSET[,"H.AFTER"]/MSET[,"AFTER"])
Spero davvero che mi possiate aiutare!!
Grazie
Laura Bettella
[[alternative HTML version deleted]]
______________________________________________ R-help at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.