[R] problemi di allocazione di memoria
Berend Hasselman
bhh at xs4all.nl
Mon Jun 17 19:06:12 CEST 2013
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.
More information about the R-help
mailing list