Luận văn Phân tích số liệu thống kê với ngôn ngữ R

models$p<-models$p[-xoa] models$d<-models$d[-xoa] models$q<-models$q[-xoa] aic<-aic[-xoa] #xep loai xl<-1:length(aic) aicxl<-sort(aic) for(k in 1:length(aic)){ for(l in 1:length(aic)){ if(aic[k]==aicxl[l]){ if(l<10){xl[k]<-paste("..........",l); break;} else if(l<100){xl[k]<-paste(".........",l); break;} }}} #tim chi so mo hinh toi uu for(t in 1:length(aic)) if(aic[t]==min(aic)) break; id.best<-t #hop nhat mo hinh mohinh<-1:length(aic) giatri.AIC<-1:length(aic) for(i in 1:length(mohinh)){ mohinh[i]<-paste("ARIMA(", models$p[i],",", models$d[i],",", models$q[i], ")",sep="") giatri.AIC[i]<-paste("...AIC =",round(aic[i],2)) } #xu ly ket qua xuat kq1<-data.frame(mohinh,giatri.AIC,xl) namescot1<-paste("Mo","hinh",sep=" ") namescot2<-paste("Gia","tri","AIC",sep=" ") namescot3<-paste("Xep","loai",sep=" ") colnames(kq1)<-c(namescot1,namescot2,namescot3) stt<-1:length(mohinh) for(j in 1:length(mohinh)) stt[j]<-paste("mo hinh",j) dimnames(kq1)[[1]]<-stt kq2<-data.frame("the.best.model"=paste(mohinh[id.best], giatri.AIC[id.best])) colnames(kq2)<-c(paste("Mo","hinh","toi","uu",sep=" ")) dimnames(kq2)[[1]]<-" " print<-list(mohinh=kq1,best=kq2) }#finish ARIMA if(type=="SARIMA"){ if (missing(order) || is.null(order) || !is.numeric(order) || length(order) != 3L || any(order < 0)) stop("'order' cau mo hinh SARIMA phai la mot vector so khong am va co chieu dai la 3!") if (is.null(seas) || !is.list(seas))stop("'seas' phai la mot danh sach gom hai thanh phan 'order' va 'frequency'!") if (is.null(seas$order) || !is.numeric(seas$order) || length(seas$order) != 3 || any(seas$order < 0)) stop("'seas$order' phai la mot vector so khong am va co chieu dai la 3!") if (is.null(seas$frequency) || is.na(seas$frequency) || seas$frequency[1] < 0 ) stop("'seas$frequency' phai la mot so nguyen khong am!") #Ghi nhan cac he so p<-order[1] d<-order[2] q<-order[3] P<-seas$order[1] D<-seas$order[2] Q<-seas$order[3] s<-seas$frequency[1] #Kiem he so mot lan nua if (is.na(p)||!is.numeric(p) || p < 0 || !is.wholenumber(p)) stop("p phai la nguyen duong") if (is.na(d)||!is.numeric(d) || d < 0 || !is.wholenumber(d)) stop("d phai la nguyen duong") if (is.na(q)||!is.numeric(q) || q < 0 || !is.wholenumber(q)) stop("p phai la nguyen duong") if (is.na(P)||!is.numeric(P) || P < 0 || !is.wholenumber(P)) stop("P phai la nguyen duong") if (is.na(D)||!is.numeric(D) || D < 0 || !is.wholenumber(D)) stop("D phai la nguyen duong") if (is.na(Q)||!is.numeric(Q) || Q < 0 || !is.wholenumber(Q)) stop("Q phai la nguyen duong") if (is.na(s)||!is.numeric(s) || s < 0 || !is.wholenumber(s)) stop("s phai la nguyen duong") #To hop mo hinh, tinh aic anser<-gl(1,1,length=(p+1)*(q+1)*(P+1)*(Q+1),labels=c("aic = ","aic =")) models<- list(p=1:((p+1)*(q+1)*(P+1)*(Q+1)),d=1:((p+1)*(q+1)*(P+1)*(Q+ 1)), q=1:((p+1)*(q+1)*(P+1)*(Q+1)),P=1:((p+1)*(q+1)*(P+1)*(Q+1)), D=1:((p+1)*(q+1)*(P+1)*(Q+1)),Q=1:((p+1)*(q+1)*(P+1)*(Q+1)), s=1:((p+1)*(q+1)*(P+1)*(Q+1))) aic<-1:((p+1)*(q+1)*(P+1)*(Q+1)) index=0 for(i1 in 0:p){ for(i2 in 0:q){ for(j1 in 0:P){ for(j2 in 0:Q){ index=index+1 models$p[index]<-i1 models$d[index]<-d models$q[index]<-i2 models$P[index]<-j1 models$D[index]<-D models$Q[index]<-j2 models$s[index]<-s aic[index]<- AIC(arima(DataTimeSeries,order=c(i1,d,i2),seasonal=list(order=c (j1,D,j2),s))) }}}} PHÂN TÍCH SỐ LIỆU THỐNG KÊ VỚI NGÔN NGỮ R 240 #Xoa bo cac mo hinh 0 co y nghia xoa<-0 for(i in 1:length(aic)) if(models$p[i]==0 & models$q[i]==0) xoa<-cbind(xoa,i) for(i in 1:length(aic)) if(models$P[i]==0 & models$Q[i]==0) xoa<-cbind(xoa,i) xoa<-unique(as.numeric(xoa[-1])) models$p<-models$p[-xoa] models$d<-models$d[-xoa] models$q<-models$q[-xoa] models$P<-models$P[-xoa] models$D<-models$D[-xoa] models$Q<-models$Q[-xoa] models$s<-models$s[-xoa] aic<-aic[-xoa] #xep loai xl<-1:length(aic) aicxl<-sort(aic) #do tim va cap so for(k in 1:length(aic)){ for(l in 1:length(aic)){ if(aic[k]==aicxl[l]){ if(l<10){xl[k]<-paste("..........",l); break;} else if(l<100){xl[k]<-paste(".........",l); break;} }}} #tim chi so mo hinh toi uu for(t in 1:length(aic)) if(aic[t]==min(aic)) break; id.best<-t #hop nhat mo hinh mohinh<-1:length(aic) giatri.AIC<-1:length(aic) for(i in 1:length(mohinh)){ mohinh[i]<- paste("SARIMA(",models$p[i],",",models$d[i],",",models$q[i], ")*(",models$P[i],",",models$D[i],",",models$Q[i], "),s=", models$s[i],sep="") giatri.AIC[i]<-paste("...AIC =",round(aic[i],2))} #xu ly ket qua in ra kq1<-data.frame(mohinh,giatri.AIC,xl) namescot1<-paste("Mo","hinh",sep=" ") namescot2<-paste("Gia","tri","AIC",sep=" ") namescot3<-paste("Xep","loai",sep=" ") colnames(kq1)<-c(namescot1,namescot2,namescot3) stt<-1:length(mohinh) for(j in 1:length(mohinh)) stt[j]<-paste("mo hinh",j) dimnames(kq1)[[1]]<-stt kq2<-data.frame("the.best.model"=paste(mohinh[id.best], giatri.AIC[id.best])) colnames(kq2)<-c(paste("Mo","hinh","toi","uu",sep=" ")) dimnames(kq2)[[1]]<-" " print<-list(mohinh=kq1,best=kq2) }#finish SARIMA if(type=="ARMAX"){ if(is.null(xreg))stop("Phai khai bao xreg cho mo hinh!") if(is.vector(xreg)) if(length(xreg)!=length(DataTimeSeries))stop("Chieu dai chuoi quan sat va xreg khac nhau!") if(is.data.frame(xreg)) if(dim(xreg)[1]!=length(DataTimeSeries))stop("Chieu dai chuoi quan sat va xreg khac nhau!" if(!is.null(seas$order)||!is.null(seas$frequency))stop("Mo hinh khong phu hop!") if (missing(order) ||is.null(order) || !is.numeric(order) || length(order) != 2L || any(order < 0)) stop("Loi o bac cua mo hinh!") #Ghi nhan cac he so p<-order[1] d<-0 q<-order[2] #Kiem he so mot lan nua if (is.na(p)||!is.numeric(p) || p < 0 || !is.wholenumber(p)) stop("Loi o bac cua mo hinh!") if (is.na(q)||!is.numeric(q) || q < 0 || !is.wholenumber(q)) stop("Loi o bac cua mo hinh!") #To hop mo hinh, tinh aic anser<-gl(1,1,length=(p+1)*(q+1),labels=c("aic = ","aic =")) models<- list(p=1:((p+1)*(q+1)),d=1:((p+1)*(q+1)),q=1:((p+1)*(q+1))) aic<-1:((p+1)*(q+1)) index=0 for(i1 in 0:p){ for(i2 in 0:q){ index=index+1 models$p[index]<-i1 models$d[index]<-d models$q[index]<-i2 aic[index]<- AIC(arimax(DataTimeSeries,order=c(i1,d,i2),xreg=xreg)) }} #Xoa bo cac mo hinh 0 co y nghia xoa<-0 for(i in 1:length(aic)) if(models$p[i]==0 & models$q[i]==0) xoa<-cbind(xoa,i) models$p<-models$p[-xoa] models$d<-models$d[-xoa] models$q<-models$q[-xoa] aic<-aic[-xoa] #xep loai xl<-1:length(aic) aicxl<-sort(aic) for(k in 1:length(aic)){ for(l in 1:length(aic)){ if(aic[k]==aicxl[l]){ if(l<10){xl[k]<-paste("..........",l); break;} else if(l<100){xl[k]<-paste(".........",l); break;} }}} #tim chi so mo hinh toi uu for(t in 1:length(aic)) if(aic[t]==min(aic)) break; id.best<-t #hop nhat mo hinh mohinh<-1:length(aic) giatri.AIC<-1:length(aic) for(i in 1:length(mohinh)){ mohinh[i]<-paste("ARMAX(", models$p[i],",", models$q[i], ")",sep="") giatri.AIC[i]<-paste("...AIC =",round(aic[i],2)) } #xu ly ket qua xuat PHÂN TÍCH SỐ LIỆU THỐNG KÊ VỚI NGÔN NGỮ R 241 kq1<-data.frame(mohinh,giatri.AIC,xl) namescot1<-paste("Mo","hinh",sep=" ") namescot2<-paste("Gia","tri","AIC",sep=" ") namescot3<-paste("Xep","loai",sep=" ") colnames(kq1)<-c(namescot1,namescot2,namescot3) stt<-1:length(mohinh) for(j in 1:length(mohinh)) stt[j]<-paste("mo hinh",j) dimnames(kq1)[[1]]<-stt kq2<-data.frame("the.best.model"=paste(mohinh[id.best], giatri.AIC[id.best])) colnames(kq2)<-c(paste("Mo","hinh","toi","uu",sep=" ")) dimnames(kq2)[[1]]<-" " print<-list(mohinh=kq1,best=kq2) }#finish ARMAX if(type=="ARIMAX"){ if(is.null(xreg))stop("Phai khai bao xreg cho mo hinh!") if(is.vector(xreg)) if(length(xreg)!=length(DataTimeSeries))stop("Chieu dai chuoi quan sat va xreg khac nhau!") if(is.data.frame(xreg)) if(dim(xreg)[1]!=length(DataTimeSeries))stop("Chieu dai chuoi quan sat va xreg khac nhau!") if(!is.null(seas$order)||!is.null(seas$frequency))stop("Mo hinh khong phu hop!") if (missing(order) ||is.null(order) || !is.numeric(order) || length(order) != 3L || any(order < 0)) stop("Loi o bac cua mo hinh!") #Ghi nhan cac he so p<-order[1] d<-order[2] q<-order[3] #Kiem he so mot lan nua if (is.na(p)||!is.numeric(p) || p < 0 || !is.wholenumber(p)) stop("Loi o bac cua mo hinh!") if (is.na(d)||!is.numeric(d) || d < 0 || !is.wholenumber(d)) stop("Loi o bac cua mo hinh!") if (is.na(q)||!is.numeric(q) || q < 0 || !is.wholenumber(q)) stop("Loi o bac cua mo hinh!") #To hop mo hinh, tinh aic anser<-gl(1,1,length=(p+1)*(q+1),labels=c("aic = ","aic =")) models<- list(p=1:((p+1)*(q+1)),d=1:((p+1)*(q+1)),q=1:((p+1)*(q+1))) aic<-1:((p+1)*(q+1)) index=0 for(i1 in 0:p){ for(i2 in 0:q){ index=index+1 models$p[index]<-i1 models$d[index]<-d models$q[index]<-i2 aic[index]<- AIC(arimax(DataTimeSeries,order=c(i1,d,i2),xreg=xreg)) }} #Xoa bo cac mo hinh 0 co y nghia xoa<-0 for(i in 1:length(aic)) if(models$p[i]==0 & models$q[i]==0) xoa<-cbind(xoa,i) models$p<-models$p[-xoa] models$d<-models$d[-xoa] models$q<-models$q[-xoa] aic<-aic[-xoa] #xep loai xl<-1:length(aic) aicxl<-sort(aic) for(k in 1:length(aic)){ for(l in 1:length(aic)){ if(aic[k]==aicxl[l]){ if(l<10){xl[k]<-paste("..........",l); break;} else if(l<100){xl[k]<-paste(".........",l); break;} }}} #tim chi so mo hinh toi uu for(t in 1:length(aic)) if(aic[t]==min(aic)) break; id.best<-t #hop nhat mo hinh mohinh<-1:length(aic) giatri.AIC<-1:length(aic) for(i in 1:length(mohinh)){ mohinh[i]<-paste("ARIMAX(", models$p[i],",", models$d[i],",", models$q[i], ")",sep="") giatri.AIC[i]<-paste("...AIC =",round(aic[i],2)) } #xu ly ket qua xuat kq1<-data.frame(mohinh,giatri.AIC,xl) namescot1<-paste("Mo","hinh",sep=" ") namescot2<-paste("Gia","tri","AIC",sep=" ") namescot3<-paste("Xep","loai",sep=" ") colnames(kq1)<-c(namescot1,namescot2,namescot3) stt<-1:length(mohinh) for(j in 1:length(mohinh)) stt[j]<-paste("mo hinh",j) dimnames(kq1)[[1]]<-stt kq2<-data.frame("the.best.model"=paste(mohinh[id.best], giatri.AIC[id.best])) colnames(kq2)<-c(paste("Mo","hinh","toi","uu",sep=" ")) dimnames(kq2)[[1]]<-" " print<-list(mohinh=kq1,best=kq2) }#finish ARIMAX if(type=="SARIMAX"){ if(is.null(xreg))stop("Phai khai bao xreg cho mo hinh!") if(is.vector(xreg)) if(length(xreg)!=length(DataTimeSeries))stop("Chieu dai chuoi quan sat va xreg khac nhau!") if(is.data.frame(xreg)) if(dim(xreg)[1]!=length(DataTimeSeries))stop("Chieu dai chuoi quan sat va xreg khac nhau!") if (missing(order) || is.null(order) || !is.numeric(order) || length(order) != 3L || any(order < 0)) stop("'order' cau mo hinh SARIMA phai la mot vector so khong am va co chieu dai la 3!") if (is.null(seas) || !is.list(seas))stop("'seas' phai la mot danh sach gom hai thanh phan 'order' va 'frequency'!") if (is.null(seas$order) || !is.numeric(seas$order) || length(seas$order) != 3 || any(seas$order < 0)) stop("'seas$order' phai la mot vector so khong am va co chieu dai la 3!") if (is.null(seas$frequency) || is.na(seas$frequency) || seas$frequency[1] < 0 ) stop("'seas$frequency' phai la mot so nguyen khong am!") PHÂN TÍCH SỐ LIỆU THỐNG KÊ VỚI NGÔN NGỮ R 242 #Ghi nhan cac he so p<-order[1] d<-order[2] q<-order[3] P<-seas$order[1] D<-seas$order[2] Q<-seas$order[3] s<-seas$frequency[1] #Kiem he so mot lan nua if (is.na(p)||!is.numeric(p) || p < 0 || !is.wholenumber(p)) stop("p phai la nguyen duong") if (is.na(d)||!is.numeric(d) || d < 0 || !is.wholenumber(d)) stop("d phai la nguyen duong") if (is.na(q)||!is.numeric(q) || q < 0 || !is.wholenumber(q)) stop("p phai la nguyen duong") if (is.na(P)||!is.numeric(P) || P < 0 || !is.wholenumber(P)) stop("P phai la nguyen duong") if (is.na(D)||!is.numeric(D) || D < 0 || !is.wholenumber(D)) stop("D phai la nguyen duong") if (is.na(Q)||!is.numeric(Q) || Q < 0 || !is.wholenumber(Q)) stop("Q phai la nguyen duong") if (is.na(s)||!is.numeric(s) || s < 0 || !is.wholenumber(s)) stop("s phai la nguyen duong") #To hop mo hinh, tinh aic anser<-gl(1,1,length=(p+1)*(q+1)*(P+1)*(Q+1),labels=c("aic = ","aic =")) models<- list(p=1:((p+1)*(q+1)*(P+1)*(Q+1)),d=1:((p+1)*(q+1)*(P+1)*(Q+ 1)), q=1:((p+1)*(q+1)*(P+1)*(Q+1)),P=1:((p+1)*(q+1)*(P+1)*(Q+1)), D=1:((p+1)*(q+1)*(P+1)*(Q+1)),Q=1:((p+1)*(q+1)*(P+1)*(Q+1)), s=1:((p+1)*(q+1)*(P+1)*(Q+1))) aic<-1:((p+1)*(q+1)*(P+1)*(Q+1)) index=0 for(i1 in 0:p){ for(i2 in 0:q){ for(j1 in 0:P){ for(j2 in 0:Q){ index=index+1 models$p[index]<-i1 models$d[index]<-d models$q[index]<-i2 models$P[index]<-j1 models$D[index]<-D models$Q[index]<-j2 models$s[index]<-s aic[index]<- AIC(arimax(DataTimeSeries,order=c(i1,d,i2),seasonal=list(order= c(j1,D,j2),s),xreg=xreg)) }}}} #Xoa bo cac mo hinh 0 co y nghia xoa<-0 for(i in 1:length(aic)) if(models$p[i]==0 & models$q[i]==0) xoa<-cbind(xoa,i) for(i in 1:length(aic)) if(models$P[i]==0 & models$Q[i]==0) xoa<-cbind(xoa,i) xoa<-unique(as.numeric(xoa[-1])) models$p<-models$p[-xoa] models$d<-models$d[-xoa] models$q<-models$q[-xoa] models$P<-models$P[-xoa] models$D<-models$D[-xoa] models$Q<-models$Q[-xoa] models$s<-models$s[-xoa] aic<-aic[-xoa] #xep loai xl<-1:length(aic) aicxl<-sort(aic) #do tim va cap so for(k in 1:length(aic)){ for(l in 1:length(aic)){ if(aic[k]==aicxl[l]){ if(l<10){xl[k]<-paste("..........",l); break;} else if(l<100){xl[k]<-paste(".........",l); break;} }}} #tim chi so mo hinh toi uu for(t in 1:length(aic)) if(aic[t]==min(aic)) break; id.best<-t #hop nhat mo hinh mohinh<-1:length(aic) giatri.AIC<-1:length(aic) for(i in 1:length(mohinh)){ mohinh[i]<- paste("SARIMAX(",models$p[i],",",models$d[i],",",models$q[i], ")*(",models$P[i],",",models$D[i],",",models$Q[i], "),s=", models$s[i],sep="") giatri.AIC[i]<-paste("...AIC =",round(aic[i],2))} #xu ly ket qua in ra kq1<-data.frame(mohinh,giatri.AIC,xl) namescot1<-paste("Mo","hinh",sep=" ") namescot2<-paste("Gia","tri","AIC",sep=" ") namescot3<-paste("Xep","loai",sep=" ") colnames(kq1)<-c(namescot1,namescot2,namescot3) stt<-1:length(mohinh) for(j in 1:length(mohinh)) stt[j]<-paste("mo hinh",j) dimnames(kq1)[[1]]<-stt kq2<-data.frame("the.best.model"=paste(mohinh[id.best], giatri.AIC[id.best])) colnames(kq2)<-c(paste("Mo","hinh","toi","uu",sep=" ")) dimnames(kq2)[[1]]<-" " print<-list(mohinh=kq1,best=kq2) }#finish SARIMAX if(type=="ARCH"){ if(!is.null(seas$order)||!is.null(seas$frequency))stop("Mo hinh khong phu hop!") if (missing(order) || is.null(order) ||!is.numeric(order) || length(order) != 1L || any(order < 0)) stop("'order' cua mo hinh ARCH phai la mot so nguyen duong") #Ghi nhan cac he so p<-0 q<-order[1] if (is.na(q)||!is.numeric(q) || q < 0 || !is.wholenumber(q)) stop("q phai la nguyen duong") #To hop mo hinh, tinh aic anser<-gl(1,1,length=(p+1)*(q+1),labels=c("aic = ","aic =")) models<-list(p=1:((p+1)*(q+1)),q=1:((p+1)*(q+1))) aic<-1:((p+1)*(q+1)) index=0 for(i1 in 0:p){ for(i2 in 0:q){ index=index+1 PHÂN TÍCH SỐ LIỆU THỐNG KÊ VỚI NGÔN NGỮ R 243 models$p[index]<-i1 models$q[index]<-i2 if(index!=1)aic[index]<- AIC(garch(DataTimeSeries,order=c(i1,i2),trace=FALSE)) }} aic<-aic[-1] models$p<-models$p[-1] models$q<-models$q[-1] #xep loai xl<-1:length(aic) aicxl<-sort(aic) for(k in 1:length(aic)){ for(l in 1:length(aic)){ if(aic[k]==aicxl[l]){ if(l<10){xl[k]<-paste("..........",l); break;} else if(l<100){xl[k]<-paste(".........",l); break;} }}} #tim chi so mo hinh toi uu for(t in 1:length(aic)) if(aic[t]==min(aic)) break; id.best<-t #Suat ket qua #hop nhat mo hinh mohinh<-1:length(aic) giatri.AIC<-1:length(aic) for(i in 1:length(aic)){ mohinh[i]<-paste("ARCH(",models$q[i], ")",sep="") giatri.AIC[i]<-paste("...AIC =",round(aic[i],2)) } #xu ly ket qua xuat kq1<-data.frame("Cac mo hinh ARCH.." = mohinh,"...Gia tri AIC.." = giatri.AIC,"....xep loai"=xl) namescot1<-paste("Mo","hinh",sep=" ") namescot2<-paste("Gia","tri","AIC",sep=" ") namescot3<-paste("Xep","loai",sep=" ") colnames(kq1)<-c(namescot1,namescot2,namescot3) stt<-1:length(mohinh) for(j in 1:length(mohinh)) stt[j]<-paste("mo hinh",j) dimnames(kq1)[[1]]<-stt kq2<-data.frame("the.best.model"=paste(mohinh[id.best], giatri.AIC[id.best])) colnames(kq2)<-c(paste("Mo","hinh","toi","uu",sep=" ")) dimnames(kq2)[[1]]<-" " print<-list(mohinh=kq1,best=kq2) }#finish ARCH if(type=="GARCH"){ if(!is.null(seas$order)||!is.null(seas$frequency))stop("Mo hinh khong phu hop!") if (missing(order) || is.null(order) ||!is.numeric(order) || length(order) != 2L || any(order < 0)) stop("'order' cua mo hinh GARCH phai la mot vector gom hai so nguyen duong") #Ghi nhan cac he so p<-order[1] q<-order[2] if (is.na(p)||!is.numeric(p) || p < 0 || !is.wholenumber(p)) stop("p phai la nguyen duong") if (is.na(q)||!is.numeric(q) || q < 0 || !is.wholenumber(q)) stop("q phai la nguyen duong") #To hop mo hinh, tinh aic anser<-gl(1,1,length=(p+1)*(q+1),labels=c("aic = ","aic =")) models<-list(p=1:((p+1)*(q+1)),q=1:((p+1)*(q+1))) aic<-1:((p+1)*(q+1)) index=0 for(i1 in 0:p){ for(i2 in 0:q){ index=index+1 models$p[index]<-i1 models$q[index]<-i2 if(i1!=0)aic[index]<- AIC(garch(DataTimeSeries,order=c(i1,i2),trace=FALSE)) }} aic<-aic[-c(1:(q+1))] models$p<-models$p[-c(1:(q+1))] models$q<-models$q[-c(1:(q+1))] #xep loai xl<-1:length(aic) aicxl<-sort(aic) for(k in 1:length(aic)){ for(l in 1:length(aic)){ if(aic[k]==aicxl[l]){ if(l<10){xl[k]<-paste("..........",l); break;} else if(l<100){xl[k]<-paste(".........",l); break;} }}} #tim chi so mo hinh toi uu for(t in 1:length(aic)) if(aic[t]==min(aic)) break; id.best<-t #Suat ket qua #hop nhat mo hinh mohinh<-1:length(aic) giatri.AIC<-1:length(aic) for(i in 1:length(aic)){ mohinh[i]<-paste("GARCH(",models$p[i],",",models$q[i], ")",sep="") giatri.AIC[i]<-paste("...AIC =",round(aic[i],2)) } #xu ly ket qua xuat kq1<-data.frame("Cac mo hinh GARCH.." = mohinh,"...Gia tri AIC.." = giatri.AIC,"....xep loai"=xl) namescot1<-paste("Mo","hinh",sep=" ") namescot2<-paste("Gia","tri","AIC",sep=" ") namescot3<-paste("Xep","loai",sep=" ") colnames(kq1)<-c(namescot1,namescot2,namescot3) stt<-1:length(mohinh) for(j in 1:length(mohinh)) stt[j]<-paste("mo hinh",j) dimnames(kq1)[[1]]<-stt kq2<-data.frame("the.best.model"=paste(mohinh[id.best], giatri.AIC[id.best])) colnames(kq2)<-c(paste("Mo","hinh","toi","uu",sep=" ")) dimnames(kq2)[[1]]<-" " print<-list(mohinh=kq1,best=kq2) }#finsh GARCH print } PHÂN TÍCH SỐ LIỆU THỐNG KÊ VỚI NGÔN NGỮ R 244 Hàm fuzzy.ts1( ) fuzzy.ts1 <- function(ts,n=5,D1=0,D2=0,type=c("Chen","Singh","Heuristic"," Chen-Hsu"),bin=NULL,trace=FALSE,divide=NULL,plot=FALSE){ #--ham con--- is.wholenumber <-function(x, tol = .Machine$double.eps^0.5) abs(x - round(x)) < tol matdo<-function(x,g){ method <- 1 x.unique <- sort(unique(x)) min.dif <- min(diff(x.unique),na.rm=TRUE)/2 min.dif.factor <- 1 nnm <- sum(!is.na(x)) n <- table(x) xx <- as.double(names(n)) cum <- cumsum(n) m <- length(xx) y <- as.integer(ifelse(is.na(x), NA, 1)) cuts <- approx(cum,xx,xout=(1:g)*nnm/g,method="constant", rule=2,f=1)$y cuts[length(cuts)] <- max(xx) lower <- xx[1] upper <- 1e+45 up <- low <- double(g) i <- 0 for (j in 1:g) { cj <- if (method == 1 || j == 1) cuts[j] else { s <- if (is.na(lower)) FALSE else xx >= lower cum.used <- if (all(s)) 0 else max(cum[!s]) if (j == m) max(xx) else if (sum(s) < 2) max(xx) else approx(cum[s] - cum.used, xx[s], xout = (nnm - cum.used)/(g - j + 1), method = "constant", rule = 2, f = 1)$y } if (cj == upper) next i <- i + 1 upper <- cj y[x >= (lower - min.dif.factor * min.dif)] <- i low[i] <- lower lower <- if (j == g) upper else min(xx[xx > upper],na.rm=TRUE) if (is.na(lower)) lower <- upper up[i] <- lower } low <- low[1:i] up <- up[1:i] kq=c(low,up[length(up)]) v=rep(0,length(kq)) for(i in 1:length(kq)) {tt=length(x[x==kq[i]]) if(tt!=0)v[i]=i} v[v==0]<-NA v<-na.omit(v) for(i in 1:length(v)) kq[v[i]]=kq[v[i]]-10^(-3) kq } rule1<-function(x2,x1,A3,D){ #trich thong tin tapmo<-as.character(D[,1]) low<-D[,2] up<-D[,3] bw<-D[,4] for(i in 1:length(tapmo)) #dinh vi tri A3 if(A3==tapmo[i]) {vt=i;break} s<-(up[vt]-low[vt])/4 #tinh khoang chia different<-abs(x2-x1)/2 distant=abs(up[vt]-low[vt])/2 if(different > distant) x=low[vt]+3*s else if(different == distant) x=low[vt]+2*s else if(different < distant) x=low[vt]+1*s x } rule2<-function(x3,x2,x1,A4,D){ #trich thong tin tapmo<-as.character(D[,1]) low<-D[,2] up<-D[,3] bw<-D[,4] for(i in 1:length(tapmo)) #dinh vi tri A2 if(A4==tapmo[i]) {vt=i;break} s<-(up[vt]-low[vt])/4 #tinh khoang chia different1a<- abs((x3-x2)-(x2-x1))*2+x3 different1b<- x3-abs((x3-x2)-(x2-x1))*2 different2a<- abs((x3-x2)-(x2-x1))/2+x3 different2b<- x3-abs((x3-x2)-(x2-x1))/2 if(low[vt]<=different1a & different1a<=up[vt]) x=low[vt]+3*s else if(low[vt]<=different1b & different1b<=up[vt]) x=low[vt]+3*s else if(low[vt]<=different2a & different2a<=up[vt]) x=low[vt]+1*s else if(low[vt]<=different2b & different2b<=up[vt]) x=low[vt]+1*s else x=low[vt]+2*s x } rule3<-function(x3,x2,x1,A4,D){ #trich thong tin tapmo<-as.character(D[,1]) low<-D[,2] up<-D[,3] bw<-D[,4] for(i in 1:length(tapmo)) #dinh vi tri A2 if(A4==tapmo[i]) {vt=i;break} s<-(up[vt]-low[vt])/4 #tinh khoang chia different1a<- abs((x3-x2)-(x2-x1))*2+x3 different1b<- x3-abs((x3-x2)-(x2-x1))*2 different2a<- abs((x3-x2)-(x2-x1))/2+x3 different2b<- x3-abs((x3-x2)-(x2-x1))/2 if(low[vt]<=different2a & different2a<=up[vt]) x=low[vt]+1*s else if(low[vt]<=different2b & different2b<=up[vt]) x=low[vt]+1*s else if(low[vt]<=different1a & different1a<=up[vt]) x=low[vt]+3*s else if(low[vt]<=different1b & different1b<=up[vt]) x=low[vt]+3*s else PHÂN TÍCH SỐ LIỆU THỐNG KÊ VỚI NGÔN NGỮ R 245 x=low[vt]+2*s x } namthang<-function(data.ts){ batdau<-start(data.ts) tanso<-frequency(data.ts) nam1<-batdau[1] thang1<-batdau[2] ketthuc<-end(data.ts) nam2<-ketthuc[1] thang2<-ketthuc[2] namkq<-1:length(data.ts) thangkq<-1:length(data.ts) index=0; for(nam in nam1:nam2) for(thang in 1:tanso) if(nam!=nam1 || thang>=thang1){ index=index+1 namkq[index]<-nam thangkq[index]<-thang if(nam==nam2 & thang==thang2) break; } if(tanso==4) {thangkq[thangkq==1]<-"Q1";thangkq[thangkq==2]<-"Q1"; thangkq[thangkq==3]<-"Q1";thangkq[thangkq==4]<-"Q1"; print<-paste(namkq,thangkq,sep=" ")} else if(tanso==12) {thangkq[thangkq==1]<-"Jan";thangkq[thangkq==2]<-"Feb"; thangkq[thangkq==3]<-"Mar";thangkq[thangkq==4]<-"Apr"; thangkq[thangkq==5]<-"May";thangkq[thangkq==6]<-"Jun"; thangkq[thangkq==7]<-"Jul";thangkq[thangkq==8]<-"Aug"; thangkq[thangkq==9]<-"Sep";thangkq[thangkq==10]<-"Oct"; thangkq[thangkq==11]<-"Nov";thangkq[thangkq==12]<-"Dec"; print<-paste(namkq,thangkq,sep=" ")} else if(tanso==7) {thangkq[thangkq==1]<-"Mon";thangkq[thangkq==2]<-"Tue"; thangkq[thangkq==3]<-"Wed";thangkq[thangkq==4]<-"Thu"; thangkq[thangkq==5]<-"Fri";thangkq[thangkq==6]<-"Sat"; thangkq[thangkq==7]<-"Sun"; print<-paste(namkq,thangkq,sep=" ")} else if(tanso!=1) print<-paste("(",namkq,",",thangkq,")",sep="") else print<-namkq print } divide.distance<-function(tsi,f,min.tsi,max.tsi){ min.x=min(tsi,min.tsi) max.x=max(tsi,max.tsi) h=(max.x-min.x)/f k<-1:(f+1); for(i in 1:(f+1)){ if(i==1)k[i]=min.x else k[i]=min.x+(i-1)*h } k } quanhe<-function(ts,type){ qh<-1:length(ts) if(type=="Chen") for(i in 1:length(ts))if(i>1)qh[i]<-c("<--") if(type=="Singh") for(i in 1:length(ts))if(i>1)qh[i]<-c("<--") if(type=="Heuristic") for(i in 1:length(ts))if(i>1)qh[i]<-c("<--") if(type=="Chen-Hsu") for(i in 1:length(ts))if(i>1)qh[i]<-c("<--") qh[qh!="<--"]<-"-x-" qh} #--ham con--- #---kiem tra input--- #Kiem tra so lieu if (!is.numeric(ts)) stop("Error in 'ts'!") if(!is.ts(ts)) stop("Error in 'ts'!") else if(!is.null(dim(ts)))stop("Error in 'ts'!") kt<-0 for(i in 1:length(ts))if(is.na(ts[i]))kt=kt+1 if(kt>0)stop("Trong chuoi co gia tri NA!") if (is.na(D1)||!is.numeric(D1)) stop("Error in 'D1'!") if (is.na(D2)||!is.numeric(D2)) stop("Error in 'D2'!") if (is.na(n)||!is.numeric(n) || n < 1 || !is.wholenumber(n)) stop("Error in 'n'!") type<-match.arg(type) if(type!="Chen" & type!="Singh" & type!="Heuristic" & type!="Chen-Hsu") stop("Error in 'type'!") if(type!="Chen-Hsu")if(!is.null(bin))stop("'bin' just for Chen-Hsu model!") if(!is.null(divide))if(type!="Chen-Hsu") stop("'divide' just for Chen-Hsu model!") if(type=="Chen-Hsu")if(is.null(divide))divide<-c("distance") else if(divide!="distance" & divide!="density")stop("'divide' must be 'distance' or 'density'!") #---kiem tra input--- #---tinh toan ban dau--- d1=D1 d2=D2 thoidiem<-namthang(ts) quanhemo<-quanhe(ts,type) #Buoc 1: Xay dung va chia tap nen min.x=min(ts)-d1 max.x=max(ts)+d2 h=(max.x-min.x)/n k<-1:(n+1);U<-1:n for(i in 1:(n+1)){ if(i==1)k[i]=min.x else {k[i]=min.x+(i-1)*h U[i-1]=paste("A",i-1,sep="")} } D<-data.frame(U,low=k[1:n],up=k[2:(n+1)]) D$Bw<-(1/2)*(D$low+D$up) #Buoc 2: Xac dinh cac tap mo loai<-1:length(ts) for(i in 1:length(ts)){ for(j in 1:n){ if (D$low[j]<=ts[i] & ts[i]<=D$up[j]) {loai[i]=paste("A",j,sep=""); break;} }} PHÂN TÍCH SỐ LIỆU THỐNG KÊ VỚI NGÔN NGỮ R 246 #Buoc 3: Xac dinh moi quan he mo loai.old<-1:length(loai) for(i in 1:length(loai)){ if(i==1)loai.old[i]=NA else loai.old[i]=loai[i-1] } quanhemokq<-paste(loai,quanhemo,loai.old,sep="") D1<-data.frame(ts,loai,loai.old) b<-table(D1$loai) ni<-1:n for(i in 1:n) for(j in 1:n){ if(paste("A",i,sep="")==names(b)[j] & j<=length(b)) {ni[i]=b[j];break;}else{ if(j>length(b)) {ni[i]=0; break;} } } D$ni<-ni #Gom nhom a<-1:(n*n) a<-matrix(a,nrow=n) for(cot in 1:n){ for(i in 2:length(D1$loai.old)){ for(j in 1:n){ if(D1$loai.old[i]==paste("A",cot,sep="") & a[j,cot]!=paste("A",j,sep="") & D1$loai[i]==paste("A",j,sep="")) {a[j,cot]=paste("A",j,sep=""); break;} }}} for(cot in 1:n){ for(j in 1:n) if(a[j,cot]!=paste("A",j,sep="")) a[j,cot]<-NA } #---tinh toan ban dau--- #---giai mo va du bao--- #type Chen if(type=="Chen"){ #1.Giai mo mo<-1:n for(cot in 1:n){ if(cot>1)a[a==1]<-NA s<-length(na.omit(a[,cot])) a[is.na(a)]<-1 if(s==1){ for(j in 1:n) if(a[j,cot]==paste("A",j,sep="")){t=D$Bw[j];break;} mo[cot]=t } else if(s>1){ t=0 for(j in 1:n) if(a[j,cot]== paste("A",j,sep=""))t=t+D$Bw[j] mo[cot]=t/s } else if(s==0){ mo[cot]=D$Bw[cot] } } #2.Du bao db<-1:length(loai.old) for(i in 1:length(loai.old)){ if(i==1)db[i]=NA else for(j in 1:n)if(loai.old[i]==paste("A",j,sep="")) {db[i]=mo[j];break;} } db<-ts(db,start=start(ts),frequency=frequency(ts)) DB1<-data.frame(thoidiem,D1$ts,quanhemokq,db) }#finish Chen #type == Singh if(type=="Singh"){ F<-1:length(ts);Ds<-1:length(ts) X<-1:length(ts);XX<-1:length(ts) Y<-1:length(ts);YY<-1:length(ts) P<-1:length(ts);PP<-1:length(ts) Q<-1:length(ts);QQ<-1:length(ts) G<-1:length(ts);GG<-1:length(ts) H<-1:length(ts);HH<-1:length(ts) E<-ts; #Xoa bo level thua: begin tt=0;v<-1:length(D$ni) for(t0 in 1:length(D$ni)) if(D$ni[t0]==0){tt=tt+1;v[tt]=t0;} if(tt==0)(v=NULL) else v=v[1:tt] if(!is.null(v)){ Dt<-D[-v,] t1=1; while(t1<=length(levels(Dt$U))){ test=0; for(t2 in 1:length(Dt$U)) if(levels(Dt$U)[t1]==Dt$U[t2]){test=1;break;} if(test==0)levels(Dt$U)[t1]=NA else t1=t1+1; }} if(is.null(v))Dt<-D #Xoa bo level thua: the end for(i in 0:(length(ts)-1)){ if(i<3)F[i+1]<-NA else{ R=0;S=0; Ds[i]=abs(abs(E[i]-E[i-1])-abs(E[i-1]-E[i-2])) X[i]=E[i]+(Ds[i]/2) XX[i]=E[i]-(Ds[i]/2) Y[i]=E[i]+Ds[i] YY[i]=E[i]-Ds[i] P[i]=E[i]+(Ds[i]/4) PP[i]=E[i]-(Ds[i]/4) Q[i]=E[i]+2*Ds[i] QQ[i]=E[i]-2*Ds[i] G[i]=E[i]+(Ds[i]/6) GG[i]=E[i]-(Ds[i]/6) H[i]=E[i]+3*Ds[i] HH[i]=E[i]-3*Ds[i] Dt.U<-as.character(Dt$U) D1.loai<-as.character(D1$loai) for(j in 1:length(Dt.U)) if(D1.loai[i+1]==Dt.U[j]) {l=j;break;} if(Dt$low[l]<=X[i] & X[i]<=Dt$up[l]) {R=R+X[i];S=S+1;} if(Dt$low[l]<=XX[i] & XX[i]<=Dt$up[l]) {R=R+XX[i];S=S+1;} PHÂN TÍCH SỐ LIỆU THỐNG KÊ VỚI NGÔN NGỮ R 247 if(Dt$low[l]<=Y[i] & Y[i]<=Dt$up[l]) {R=R+Y[i];S=S+1;} if(Dt$low[l]<=YY[i] & YY[i]<=Dt$up[l]) {R=R+YY[i];S=S+1;} if(Dt$low[l]<=P[i] & P[i]<=Dt$up[l]) {R=R+P[i];S=S+1;} if(Dt$low[l]<=PP[i] & PP[i]<=Dt$up[l]) {R=R+PP[i];S=S+1;} if(Dt$low[l]<=Q[i] & Q[i]<=Dt$up[l]) {R=R+Q[i];S=S+1;} if(Dt$low[l]<=QQ[i] & QQ[i]<=Dt$up[l]) {R=R+QQ[i];S=S+1;} if(Dt$low[l]<=G[i] & G[i]<=Dt$up[l]) {R=R+G[i];S=S+1;} if(Dt$low[l]<=GG[i] & GG[i]<=Dt$up[l]) {R=R+GG[i];S=S+1;} if(Dt$low[l]<=H[i] & H[i]<=Dt$up[l]) {R=R+H[i];S=S+1;} if(Dt$low[l]<=HH[i] & HH[i]<=Dt$up[l]) {R=R+HH[i];S=S+1;} F[i+1]=(R + Dt$Bw[l])/(S+1) } } F<-ts(F,start=start(ts),frequency=frequency(ts)) DB2<- data.frame(thoidiem,sl.goc=D1$ts,quanhemo=quanhemokq,..... sl.mo=F) }#finish Singh #type Heuristic if(type=="Heuristic"){ #Xac dinh nhom quan he mo Heuristic tap="begin";t.thai="calm";Bw.h=0 for(cot in 1:n){ locate<-rep(0,n) for(i in 1:n) if(!is.na(a[i,cot]))locate[i]=i locate[locate==0]<-NA locate<-na.omit(locate) if(length(locate)>0){ #xu huong giam: lon->nho l=0 locate[length(locate)+1]<-(n+1) for(i in 1:(length(locate)-1)) if(cot >= locate[i] & cot < locate[i+1]) {l=i;break;} locate=locate[1:(length(locate)-1)] if(l!=0){ t=0 for(i in 1:l) t=t+D$Bw[locate[i]] tb<-t/length(locate[1:l]) tap[length(tap)+1]=paste("A",cot,sep="") t.thai[length(t.thai)+1]="down" Bw.h[length(Bw.h)+1]=tb } #xu huong tang: nho->lon l=0 for(i in 1:length(locate)) if(cot <= locate[i]) {l=i;break;} if(l!=0){ t=0 for(i in l:length(locate)) t=t+D$Bw[locate[i]] tb<-t/length(locate[l:length(locate)]) tap[length(tap)+1]=paste("A",cot,sep="") t.thai[length(t.thai)+1]="raise" Bw.h[length(Bw.h)+1]=tb }}} tap<-factor(tap[2:length(tap)]) t.thai<-t.thai[2:length(t.thai)] Bw.h<-Bw.h[2:length(Bw.h)] D.h<-data.frame(tap,t.thai,Bw.h) #Du bao Heuristic<-1:length(ts) Heuristic[1]<-NA for(j in 2:length(ts)){ xj<-ts[j]-ts[j-1] D1.loai.old<-as.character(D1$loai.old) D.h.tap<-as.character(D.h$tap) if(xj>=0){ for(v in 1:dim(D.h)[1]) if(D1.loai.old[j]==D.h.tap[v] & D.h$t.thai[v]=="raise") bt=v Heuristic[j]=D.h$Bw.h[bt] } else if(xj<0){ for(v in 1:dim(D.h)[1]) if(D1.loai.old[j]==D.h.tap[v] & D.h$t.thai[v]=="down") bt=v Heuristic[j]=D.h$Bw.h[bt] } } Heuristic<-ts(Heuristic,start=start(ts),frequency=frequency(ts)) DB3<- data.frame(thoidiem,sl.goc=D1$ts,quanhemo=quanhemokq,..... sl.mo=Heuristic) }#finish Heuristic #type Chen va Hsu if(type=="Chen-Hsu"){ if(is.null(bin)){ #bin==null #Phan lai tap mo b<-D$ni level=max(round(length(ts)/n-0.5),min(b),2)#tuy theo y nguoi dung f<-1:n for(i in 1:n) {if(b[i]>(level+2))f[i]=round(b[i]/level+0.5) else f[i]=1; if(f[i]==0) f[i]=1} n2=sum(f) U2<-1:n2; k2<-1:(n2+1); t=0; k2[1]<-k[1] for( i in 1:n){ tsi=ts[ts>=D$low[i]] tsi=tsi[tsi<=D$up[i]] if(divide=="density")new=matdo(tsi,f[i]) else if(divide=="distance")new=divide.distance(tsi,f[i],D$low[i],D$up [i]) else stop("divide phai co gia tri la 'distance' hoac 'frequency'!") new[1]=min(new[1],D$low[i]) new[length(new)]=max(new[length(new)],D$up[i]) for(j in 2:(f[i]+1)){ t=t+1 k2[t+1]=new[j] U2[t]=paste("A",t,sep="") } } } else{ #bin !=null if(min(bin)=max.x || length(bin)<1) stop("gia tri tam bin khong hop le!") k2<-1:(length(bin)+2);U2<-1:(length(bin)+1) k2[1]<-min.x;k2[length(k2)]<-max.x for(t in 2:(length(k2)-1)) k2[t]=bin[t-1] for(t in 1:(length(bin)+1)) U2[t]=paste("A",t,sep="") n2<-length(U2) } PHÂN TÍCH SỐ LIỆU THỐNG KÊ VỚI NGÔN NGỮ R 248 D.ch<-data.frame(U2,low2=k2[1:n2],up2=k2[2:(n2+1)]) D.ch$Bw2<-(1/2)*(D.ch$low+D.ch$up) #Xac dinh moi quan he mo moi loai2<-1:length(ts) for(i in 1:length(ts)){ for(j in 1:n2){ if (D.ch$low2[j]<=ts[i] & ts[i]<=D.ch$up2[j]) {loai2[i]=paste("A",j,sep=""); break;} }} D1.ch<-data.frame(ts,loai2,loai2.old=c(NA,loai2[1:(length(loai2)- 1)])) b<-table(D1.ch$loai2) quanhemokq<- paste(D1.ch$loai2,quanhemo,D1.ch$loai2.old,sep="") ni<-1:n2 for(i in 1:n2) for(j in 1:n2){ if(paste("A",i,sep="")==names(b)[j] & j<=length(b)) {ni[i]=b[j]; break;}else{ if(j>length(b)){ni[i]=0;break;} } } D.ch$ni<-ni #Xoa bo level thua: begin tt=0;v<-1:length(D.ch$ni) for(t0 in 1:length(D.ch$ni)) if(D.ch$ni[t0]==0){tt=tt+1;v[tt]=t0;} if(tt==0)(v=NULL) else v=v[1:tt] if(!is.null(v)){ Dt<-D.ch[-v,] t1=1; while(t1<=length(levels(Dt$U2))){ test=0; for(t2 in 1:length(Dt$U2)) if(levels(Dt$U2)[t1]==Dt$U2[t2]){test=1;break;} if(test==0)levels(Dt$U2)[t1]=NA else t1=t1+1; }} else Dt<-D.ch #Xoa bo level thua: the end #Giai mo va du bao X=ts A<-as.character(D1.ch$loai2) P<-1:length(X) tapmochonam2<-as.character(Dt[,1]) for(i in 1:length(tapmochonam2)) #dinh vi tri tap mo cua nam 2 if(A[2]==tapmochonam2[i]) {vt=i;break} P[1]<-NA;P[2]=Dt[vt,4] for(t in 2:(length(ts)-1)){ if(t==2) P[t+1]=rule1(X[t],X[t-1],A[t+1],Dt) else { different=(X[t]-X[t-1])-(X[t-1]-X[t-2]) if(A[t+1]>A[t] & different>=0) P[t+1]=rule2(X[t],X[t-1],X[t- 2],A[t+1],Dt)#th1 else if(A[t+1]>A[t] & different<=0) P[t+1]=rule3(X[t],X[t-1],X[t- 2],A[t+1],Dt)#th2 else if(A[t+1]=0) P[t+1]=rule2(X[t],X[t-1],X[t- 2],A[t+1],Dt)#th3 else if(A[t+1]<A[t] & different<=0) P[t+1]=rule3(X[t],X[t-1],X[t- 2],A[t+1],Dt)#th4 else if(A[t+1]==A[t] & different>=0) P[t+1]=rule2(X[t],X[t-1],X[t- 2],A[t+1],Dt)#th5 else if(A[t+1]==A[t] & different<=0) P[t+1]=rule3(X[t],X[t-1],X[t- 2],A[t+1],Dt)#th6 } } P<-ts(P,start=start(ts),frequency=frequency(ts)) DB4<- data.frame(thoidiem,sl.goc=D1.ch$ts,quanhemo=quanhemokq, ".....sl.mo"=P) diem25<-1:length(D.ch[,3]) diem75<-1:length(D.ch[,3]) h<-(D.ch[,3]-D.ch[,2])/4 canduoi<-D.ch[,2] for(i in 1:length(D.ch[,3])){ diem25[i]<-canduoi[i]+1*h[i] diem75[i]<-canduoi[i]+3*h[i]} D.ch<-data.frame(D.ch[,1],D.ch[,2],diem25,D.ch[,4],diem75, D.ch[,3],D.ch[,5]) namescot1<-"set" namescot2<-"dow" namescot3<-"0.25" namescot4<-"0.50" namescot5<-"0.75" namescot6<-"up" namescot7<-"num" colnames(D.ch)<- c(namescot1,namescot2,namescot3,namescot4,namescot5, namescot6,namescot7) }#finish chen-hsu #---giai mo va du bao--- #---tinh toan do chinh xac--- if(type=="Chen") accuracy<- av.res(Y=data.frame(DB1[,2]),F=data.frame(Chen=DB1[,4])) if(type=="Singh") accuracy<- av.res(Y=data.frame(DB2[,2]),F=data.frame(Singh=DB2[,4])) if(type=="Heuristic") accuracy<- av.res(Y=data.frame(DB3[,2]),F=data.frame(Heuristic=DB3[,4])) if(type=="Chen-Hsu") accuracy<- av.res(Y=data.frame(DB4[,2]),F=data.frame("Chen- Hsu"=DB4[,4])) #---tinh toan do chinh xac--- #---xuat ket qua--- if(trace==TRUE){ namescot1<-"set" namescot2<-"dow" namescot3<-"up" namescot4<-"mid" namescot5<-"num" colnames(D)<- c(namescot1,namescot2,namescot3,namescot4,namescot5) namescot1<-"point" PHÂN TÍCH SỐ LIỆU THỐNG KÊ VỚI NGÔN NGỮ R 249 namescot2<-"ts" namescot3<-"relative" namescot4<-"forecast" if(type=="Chen")colnames(DB1)<- c(namescot1,namescot2,namescot3,namescot4) if(type=="Singh")colnames(DB2)<- c(namescot1,namescot2,namescot3,namescot4) if(type=="Heuristic")colnames(DB3)<- c(namescot1,namescot2,namescot3,namescot4) if(type=="Chen-Hsu")colnames(DB4)<- c(namescot1,namescot2,namescot3,namescot4) if(type=="Chen") MO<-list(type="Chen",table1=D,table2=DB1, accuracy=accuracy) if(type=="Singh") MO<-list(type="Singh",table1=D,table2=DB2, accuracy=accuracy) if(type=="Heuristic") MO<-list(type="Heuristic",table1=D, table2=DB3,accuracy=accuracy) if(type=="Chen-Hsu") { if(is.null(bin)) MO<-list(type="Chen-Hsu",table1=D,table2=D.ch, table3=DB4,accuracy=accuracy) else if(!is.null(bin)) MO<-list(type="Chen-Hsu",table1=D.ch, table2=DB4,accuracy=accuracy) }} else if(trace==FALSE){ if(type=="Chen") MO<-DB1[,4] if(type=="Singh") MO<-DB2[,4] if(type=="Heuristic") MO<-DB3[,4] if(type=="Chen-Hsu") MO<-DB4[,4]} else MO<-c("trace phai co gia tri 'TRUE' hoac 'FALSE'") #---xuat ket qua--- #---ve bieu do du bao--- if(plot==TRUE){ if(type=="Chen") { goc<-DB1[,2] dubao<-DB1[,4] if(length(goc)<50){ plot(goc,col="blue",main=paste("Chen",n,"fuzzy set"),type="o", pch=16,ylim=c(min(c(goc,dubao),na.rm=1),max(c(goc,dubao), na.rm=1)),xlab="index",ylab="data"); lines(dubao,col="red",type="o",pch=18) legend("bottomright","(x,y)",c("ts","forecast"),col=c("blue","red "),lty=c(1,1),pch=c(16,18)) } if(length(goc)>49){ plot(goc,col="blue",main=paste("Chen",n,"fuzzy set"),type="l", ylim=c(min(c(goc,dubao),na.rm=1),max(c(goc,dubao),na.rm=1)), xlab="index",ylab="data"); lines(dubao,col="red",type="l") legend("bottomright","(x,y)",c("ts","forecast"),col=c("blue","red "),lty=c(1,1)) } k.ve<-par()$yaxp h=(k.ve[2]-k.ve[1])/k.ve[3] for(i in 0:k.ve[3]) abline(h=k.ve[1]+h*i,lty=3,col="gray") } if(type=="Singh"){ goc<-DB2[,2] dubao<-DB2[,4] if(length(goc)<50){ plot(goc,col="blue",main=paste("Singh",n,"fuzzy set"),type="o", pch=16,ylim=c(min(c(goc,dubao),na.rm=1),max(c(goc,dubao), na.rm=1)),xlab="index",ylab="data"); lines(dubao,col="red",type="o",pch=18) legend("bottomright","(x,y)",c("ts","forecast"),col=c("blue", "red"),lty=c(1,1),pch=c(16,18)) } if(length(goc)>49){ plot(goc,col="blue",main=paste("Singh",n,"fuzzy set"),type="l", ylim=c(min(c(goc,dubao),na.rm=1),max(c(goc,dubao),na.rm=1)), xlab="index",ylab="data"); lines(dubao,col="red",type="l") legend("bottomright","(x,y)",c("ts","forecast"),col=c("blue","red "),lty=c(1,1)) } k.ve<-par()$yaxp h=(k.ve[2]-k.ve[1])/k.ve[3] for(i in 0:k.ve[3]) abline(h=k.ve[1]+h*i,lty=3,col="gray") } if(type=="Heuristic") { goc<-DB3[,2] dubao<-DB3[,4] if(length(goc)<50){ plot(goc,col="blue",main=paste("Heuristic",n,"fuzzy set"),type="o",pch=16,ylim=c(min(c(goc,dubao),na.rm=1), max(c(goc,dubao),na.rm=1)),xlab="index",ylab="data"); lines(dubao,col="red",type="o",pch=18) legend("bottomright","(x,y)",c("ts","forecast"),col=c("blue","red "),lty=c(1,1),pch=c(16,18)) } if(length(goc)>49){ plot(goc,col="blue",main=paste("Heuristic",n,"fuzzy set"), type="l",ylim=c(min(c(goc,dubao),na.rm=1),max(c(goc,dubao), na.rm=1)),xlab="index",ylab="data"); lines(dubao,col="red",type="l") legend("bottomright","(x,y)",c("ts","forecast"),col=c("blue", "red"),lty=c(1,1)) } k.ve<-par()$yaxp h=(k.ve[2]-k.ve[1])/k.ve[3] for(i in 0:k.ve[3]) abline(h=k.ve[1]+h*i,lty=3,col="gray") } if(type=="Chen-Hsu") { goc<-DB4[,2] dubao<-DB4[,4] if(length(goc)<50){ plot(goc,col="blue",main=paste("Chen-Hsu",n, "fuzzy set"),type="o", pch=16,ylim=c(min(c(goc,dubao), na.rm=1),max(c(goc,dubao),na.rm=1)),xlab="index", ylab="data"); lines(dubao,col="red",type="o",pch=18) legend("bottomright","(x,y)",c("ts","forecast"),col=c("blue", "red"),lty=c(1,1),pch=c(16,18)) } if(length(goc)>49){ PHÂN TÍCH SỐ LIỆU THỐNG KÊ VỚI NGÔN NGỮ R 250 plot(goc,col="blue",main=paste("Chen-Hsu",n,"fuzzy set"),type="l",ylim=c(min(c(goc,dubao),na.rm=1),max(c(goc, dubao),na.rm=1)),xlab="index",ylab="data"); lines(dubao,col="red",type="l") legend("bottomright","(x,y)",c("ts","forecast"),col=c("blue", "red"),lty=c(1,1)) } k.ve<-par()$yaxp h=(k.ve[2]-k.ve[1])/k.ve[3] for(i in -2:(k.ve[3]+2)) abline(h=k.ve[1]+h*i,lty=3,col="gray") } } #---ve bieu do du bao--- MO } Hàm fuzzy.ts2( ) fuzzy.ts2 <-function(ts,n=5,w=NULL,D1=0,D2=0,C=NULL,r=4, trace=FALSE,forecast=NULL,plot=FALSE){ #---ham con--- is.wholenumber <-function(x, tol = .Machine$double.eps^0.5) abs(x - round(x)) < tol namthang<-function(data.ts){ batdau<-start(data.ts) tanso<-frequency(data.ts) nam1<-batdau[1] thang1<-batdau[2] ketthuc<-end(data.ts) nam2<-ketthuc[1] thang2<-ketthuc[2] namkq<-1:length(data.ts) thangkq<-1:length(data.ts) index=0; for(nam in nam1:nam2) for(thang in 1:tanso) if(nam!=nam1 || thang>=thang1){ index=index+1 namkq[index]<-nam thangkq[index]<-thang if(nam==nam2 & thang==thang2) break; } if(tanso==4) {thangkq[thangkq==1]<-"Q1";thangkq[thangkq==2]<-"Q1"; thangkq[thangkq==3]<-"Q1";thangkq[thangkq==4]<-"Q1"; print<-paste(namkq,thangkq,sep=" ")} else if(tanso==12) {thangkq[thangkq==1]<-"Jan";thangkq[thangkq==2]<-"Feb"; thangkq[thangkq==3]<-"Mar";thangkq[thangkq==4]<-"Apr"; thangkq[thangkq==5]<-"May";thangkq[thangkq==6]<-"Jun"; thangkq[thangkq==7]<-"Jul";thangkq[thangkq==8]<-"Aug"; thangkq[thangkq==9]<-"Sep";thangkq[thangkq==10]<-"Oct"; thangkq[thangkq==11]<-"Nov";thangkq[thangkq==12]<- "Dec"; print<-paste(namkq,thangkq,sep=" ")} else if(tanso==7) {thangkq[thangkq==1]<-"Mon";thangkq[thangkq==2]<- "Tue"; thangkq[thangkq==3]<-"Wed";thangkq[thangkq==4]<-"Thu"; thangkq[thangkq==5]<-"Fri";thangkq[thangkq==6]<-"Sat"; thangkq[thangkq==7]<-"Sun"; print<-paste(namkq,thangkq,sep=" ")} else if(tanso!=1) print<-paste("(",namkq,",",thangkq,")",sep="") else print<-namkq print } computeVi<-function(M,table,w){ n<-(dim(M)[1] - w) cot<-dim(M)[2] Vi<-1:dim(M)[1] Vi[1:w]<-NA for(i in 1:n){ O<-M[i:(w-1+(i-1)),] if(w==2)O<-t(as.matrix(O)) K<-M[w+(i-1),] R<-O for (i1 in 1:cot) for (j1 in 1:(w-1)) if (O[j1,i1] > K[i1]) R[j1,i1] <- K[i1] F <- 1:cot for (i2 in 1:cot) F[i2] <- max(R[, i2]) Vi[w+i]<- sum(F * table$Bw)/sum(F) } Vi } computeVi2<-function(M,table){ cot<-dim(M)[2] dong<-(dim(M)[1]-1) O<-M[1:dong,] if(w==2)O<-t(as.matrix(O)) K<-M[(dong+1),] R<-O for (i1 in 1:cot) for (j1 in 1:(w-1)) if (O[j1,i1] > K[i1]) R[j1,i1] <- K[i1] F <- 1:cot for (i2 in 1:cot) F[i2] <- max(R[, i2]) Vi<- sum(F * table$Bw)/sum(F) Vi } #---ham con--- #---kiem tra input--- if (!is.numeric(ts)) stop("Error in 'ts'!") if(!is.ts(ts)) stop("Error in 'ts'!") else if(!is.null(dim(ts)))stop("Error in 'ts'!") kt<-0 for(i in 1:length(ts))if(is.na(ts[i]))kt=kt+1 PHÂN TÍCH SỐ LIỆU THỐNG KÊ VỚI NGÔN NGỮ R 251 if(kt>0)stop("Time series contain NA!") if (is.na(n)||!is.numeric(n) || n < 1 || !is.wholenumber(n)) stop("Error in 'n'!") if (is.na(r)||!is.numeric(r) || r < 0 || !is.wholenumber(r)) stop("Error in 'r'!") if (is.null(w)||is.na(w)||!is.numeric(w) || w < 2 || !is.wholenumber(w)) stop("Error in 'w'!") if (is.null(C)||is.na(C)||!is.numeric(C)) stop("Error in 'C'!") if (is.na(D1)||!is.numeric(D1)) stop("Error in 'D1'!") if (is.na(D2)||!is.numeric(D2)) stop("Error in 'D2'!") if (is.null(forecast)||is.na(forecast)||!is.numeric(forecast) || forecast < 1 || !is.wholenumber(forecast)) stop("Error in 'forecast'!") if(w>=length(ts)) stop("Error in 'w'!") #---kiem tra input--- #---tap bien doi--- ts1<-as.vector(diff(ts)) min.x=min(ts1)-D1 max.x=max(ts1)+D2 h=(max.x-min.x)/n k<-1:(n+1);U<-1:n for(i in 1:(n+1)){ if(i==1)k[i]=min.x else {k[i]=min.x+(i-1)*h U[i-1]=paste("u",i-1,sep="")}} D<-data.frame(U,low=k[1:n],up=k[2:(n+1)]) D$Bw<-(1/2)*(D$low+D$up) table1<-D #---tap bien doi--- #---mo hoa bien doi--- thoidiem <- namthang(ts) Ai <- 1:length(ts1) MATRIX <- matrix(1:(length(ts1) * n), ncol = n) for (i in 1:length(ts1)) { temp = "" for (j in 1:n) { my.At <- 1/(1 + (C * (ts1[i] - D$Bw[j]))^2) MATRIX[i, j] <- my.At At.j <- paste("(",round(my.At,r),"/u",j,sep = "", ")") if (j == 1) temp <- paste(temp, At.j, sep = "") else temp <- paste(temp, At.j, sep = ",") } Ai[i] <- paste("A[",thoidiem[i+1],"]={",temp,"}",sep="") } table2<-data.frame(point=thoidiem,ts=ts,diff.ts=c(NA,ts1)) table3<-c(NA,Ai) #---mo hoa bien doi--- #---noi suy--- V<-computeVi(MATRIX,table1,w) N<-c(NA,(ts[-length(ts)]+V)) V<-c(NA,V) table4<- data.frame(point=thoidiem,interpolate=N,diff.interpolate=V ) accuracy<- av.res(Y=data.frame(ts),F=data.frame("Abbasov.Mamedova" =table4[,2])) table4<-na.omit(table4) rownames(table4)<-c(1:dim(table4)[1]) #---noi suy--- #---du bao--- V<-1:forecast N<-0:forecast Ai<-1:forecast N[1]<-table4[dim(table4)[1],2] MT<-MATRIX[(dim(MATRIX)[1]-w+1):dim(MATRIX)[1],] temp<-c(as.vector(ts),1:forecast) temp<-ts(temp,start=start(ts),frequency=frequency(ts)) temp<-namthang(temp) temp<-temp[(length(ts)+1):length(temp)] thoidiem<-temp for(i in 1:forecast){ V[i]<-computeVi2(MT,table1) N[i+1]<-N[i]+V[i] for(chuyen in 1:(dim(MT)[1]-1)) MT[chuyen,]<- MT[(chuyen+1),] temp = "" for (j in 1:n) { my.At <- 1/(1 + (C * (V[i] - D$Bw[j]))^2) MT[dim(MT)[1], j] <- my.At At.j <- paste("(",round(my.At,r),"/u",j,sep = "", ")") if (j == 1) temp <- paste(temp, At.j, sep = "") else temp <- paste(temp, At.j, sep = ",") } Ai[i] <- paste("A[",thoidiem[i],"]={",temp,"}",sep="") } N<-N[-1] table5<- data.frame(point=thoidiem,forecast=N,diff.forecast=V) table6<-Ai #---du bao--- KQ<-list(type="Abbasov- Manedova",table1=table1,table2=table2,table3=table3,table 4=table4,table5=table5,table6=table6,accuracy=accuracy) if(trace==TRUE) MO<-KQ else if(trace==FALSE) MO<-table5 else MO<-c("'trace' must be 'TRUE' or 'FALSE'") if(plot==TRUE) { tsp<-ts tsp[1:w]<-NA tsp<-na.omit(tsp) goc<-ts(table2[,2],start=start(ts),frequency=frequency(ts)) dubao<- ts(c(table4[,2],table5[,2]),start=start(tsp),frequency=frequen cy(tsp)) if(length(goc)<50){ plot(goc,col="blue",main=paste("Abbasov-Mamedova, C =",C,"w =",w,"n =",n),type="o",pch=16,xlab="index", ylab="data"); lines(dubao,col="red",type="o",pch=18) PHÂN TÍCH SỐ LIỆU THỐNG KÊ VỚI NGÔN NGỮ R 252 legend("bottomright","(x,y)",c("ts","forecast"),col=c("blue"," red"),lty=c(1,1),pch=c(16,18)) } if(length(goc)>49){ plot(goc,col="blue",main=paste("Abbasov-Mamedova, C =",C,"w =",w,"n =",n),type="l",xlab="index",ylab="data"); lines(dubao,col="red",type="l") legend("bottomright","(x,y)",c("ts","forecast"),col=c("blue"," red"),lty=c(1,1)) } k.ve<-par()$yaxp h=(k.ve[2]-k.ve[1])/k.ve[3] for(i in -2:(k.ve[3]+2)) abline(h=k.ve[1]+h*i,lty=3,col="gray") } MO }