Диссертация (1151123), страница 20
Текст из файла (страница 20)
Реализация алгоритма оцениванияхарактера нескольких внешних событий сиспользованием процедуры auto-arima в пакете Rrm(list=ls(all=TRUE))data <- read.table("yota_20130114.csv", header=T, sep=";",dec=",")yota <- data$revenueyota_ts <- ts(yota, frequency=12, start=c(2009,6))#data#logcars<-log(carsts)#plot.ts(logcars)library(forecast)cwp <- function (object){## cwp <--> ``coefficients with p-values''#coef <- coef(object)if (length(coef) > 0) {mask <- object$masksdev <- sqrt(diag(vcov(object)))t.rat <- rep(NA, length(mask))t.rat[mask] <- coef[mask]/sdevpt <- 2 * pnorm(-abs(t.rat))setmp <- rep(NA, length(mask))setmp[mask] <- sdevsum <- rbind(coef, setmp, t.rat, pt)dimnames(sum) <- list(c("coef", "s.e.", "t ratio", "pvalue"),names(coef))return(sum)} else return(NA)}################################################################adf.test(yota_ts)#CASE (b,b)##1.
Simple ARIMAres.1<-arima(yota_ts, xreg=data[,3:4], order = c(1,0,0),method="ML")#acf(res.1$residuals, lag.max=20)# ACF plot#pacf(res.1$residuals, lag.max=20) # PACF plot#sum(res.1$residuals*res.1$residuals)/length(res.1$residuals)#res.1fin<c('[b,b]','simple',res.1$arma,sum(res.1$residuals*res.1$residuals)/length(res.1$residuals),res.1$coef[3],res.1$coef[4])names(fin)<c('intervention_type','model','AR','MA','SAR','SMA','period','difNonSeas','difSeas','MSE','tarifs','swap')139##2. Auto-arimares.1<-auto.arima(yota_ts, xreg=data[,3:4],max.p = 4, max.q = 4,max.P = 0, max.Q = 0)#sum(res.1$residuals*res.1$residuals)/length(res.1$residuals)#res.1fin<rbind(fin,c('[b,b]','auto',res.1$arma,sum(res.1$residuals*res.1$residuals)/length(res.1$residuals),res.1$coef[2],res.1$coef[3]))##3.
Complex-arimares.2<-arima(yota_ts, xreg=data[,3:4], order = c(1,0,1),method="ML")#acf(res.1$residuals, lag.max=20)# ACF plot#pacf(res.1$residuals, lag.max=20) # PACF plot#sum(res.2$residuals*res.2$residuals)/length(res.2$residuals)#res.2cwp(res.2)fin<rbind(fin,c('[b,b]','complex',res.2$arma,sum(res.2$residuals*res.2$residuals)/length(res.2$residuals),res.2$coef[4],res.2$coef[5]))#################################################################CASE (b,c)##1. Simple ARIMA#Case BC ready for (1,0,0)(1,0,0)a1<-0.8interv<-data.frame(row(as.matrix(data[,4])),data[,3:4])names(interv)<-c('n','var1','var2')interv$bc<-interv$var2interv$bc[interv$n>min(interv$n[interv$var2>0])]<-1a1*(interv$n[interv$n>min(interv$n[interv$var2>0])]min(interv$n[interv$var2>0]))/(max(interv$n)min(interv$n[interv$var2>0])+1)res.1<-arima(yota_ts, xreg=data.frame(interv$var1,interv$bc),order = c(1,0,0), method="ML")#res.1#sum(res.1$residuals*res.1$residuals)/length(res.1$residuals)fin<rbind(fin,c('[b,c]','simple',res.1$arma,sum(res.1$residuals*res.1$residuals)/length(res.1$residuals),res.1$coef[3],res.1$coef[4]))##2.
Auto-arimainterv<-data.frame(row(as.matrix(data[,4])),data[,3:4])names(interv)<-c('n','var1','var2')a1<-seq(0.1, 1, 0.05)zmin<-1000000000140f<-0for(i in a1) {interv$bc<-interv$var2interv$bc[interv$n>min(interv$n[interv$var2>0])]<-1i*(interv$n[interv$n>min(interv$n[interv$var2>0])]min(interv$n[interv$var2>0]))/(max(interv$n)min(interv$n[interv$var2>0])+1)res.1<-auto.arima(yota_ts,xreg=data.frame(interv$var1,interv$bc),max.p = 4, max.q = 4,max.P = 0, max.Q = 0)z<-sum(res.1$residuals*res.1$residuals)/length(res.1$residuals)if(z<zmin){zmin<-za_out<-imodelmin<-res.1interv_bc<-data.frame(interv$var1,interv$bc)}f<-f+1}fin<rbind(fin,c('[b,c]','auto',modelmin$arma,zmin,modelmin$coef[2],modelmin$coef[3]))zmina_outmodelmin#fwrite.table(interv_bc,file="interv_bc.csv",sep=";",row.names=F,dec=",")##3.
Complex-arimainterv<-data.frame(row(as.matrix(data[,4])),data[,3:4])a1<-0.8names(interv)<-c('n','var1','var2')interv$bc<-interv$var2interv$bc[interv$n>min(interv$n[interv$var2>0])]<-1a1*(interv$n[interv$n>min(interv$n[interv$var2>0])]min(interv$n[interv$var2>0]))/(max(interv$n)min(interv$n[interv$var2>0])+1)res.1<-arima(yota_ts, xreg=data.frame(interv$var1,interv$bc),order = c(1,0,1), method="ML")#acf(res.1$residuals, lag.max=20)# ACF plot#pacf(res.1$residuals, lag.max=20) # PACF plot#sum(res.1$residuals*res.1$residuals)/length(res.1$residuals)#res.1cwp(res.1)fin<rbind(fin,c('[b,c]','complex',res.1$arma,sum(res.1$residuals*res.1$residuals)/length(res.1$residuals),res.1$coef[4],res.1$coef[5]))################################################################141#CASE (c,b)##1. Simple ARIMA#Case CB ready for (1,0,0)(1,0,0)a2<-1interv<-data.frame(row(as.matrix(data[,4])),data[,3:4])names(interv)<-c('n','var1','var2')interv$cb<-interv$var1interv$cb[interv$n>min(interv$n[interv$var1>0])]<-1a2*(interv$n[interv$n>min(interv$n[interv$var1>0])]min(interv$n[interv$var1>0]))/(max(interv$n)min(interv$n[interv$var1>0])+1)res.1<-arima(yota_ts, xreg=data.frame(interv$cb,interv$var2),order = c(1,0,0), method="ML")#res.1#sum(res.1$residuals*res.1$residuals)/length(res.1$residuals)fin<rbind(fin,c('[c,b]','simple',res.1$arma,sum(res.1$residuals*res.1$residuals)/length(res.1$residuals),res.1$coef[3],res.1$coef[4]))##2.
Auto-arimainterv<-data.frame(row(as.matrix(data[,4])),data[,3:4])names(interv)<-c('n','var1','var2')zmin<-1000000000a1<-seq(0.1, 1, 0.05)for(i in a1) {interv$cb<-interv$var1interv$cb[interv$n>min(interv$n[interv$var1>0])]<-1i*(interv$n[interv$n>min(interv$n[interv$var1>0])]min(interv$n[interv$var1>0]))/(max(interv$n)min(interv$n[interv$var1>0])+1)res.1<-auto.arima(yota_ts,xreg=data.frame(interv$cb,interv$var2),max.p = 4, max.q = 4,max.P = 0, max.Q = 0)z<-sum(res.1$residuals*res.1$residuals)/length(res.1$residuals)if(z<zmin){zmin<-za_out<-imodelmin<-res.1interv_cb<-data.frame(interv$cb,interv$var2)}}zmina_outmodelminwrite.table(interv_cb,file="interv_cb.csv",sep=";",row.names=F,dec=",")fin<rbind(fin,c('[c,b]','auto',modelmin$arma,zmin,modelmin$coef[2],modelmin$coef[3]))142#interv_cb_2<-data.frame(interv_cb[,1]*(27.8092),interv_cb[,2]*(-30.86))#plot(yota_ts-interv_cb_2[,1]-interv_cb_2[,2])#lines(yota_ts, col='green')##3.
Complex-arimaa2<-1interv<-data.frame(row(as.matrix(data[,4])),data[,3:4])names(interv)<-c('n','var1','var2')interv$cb<-interv$var1interv$cb[interv$n>min(interv$n[interv$var1>0])]<-1a2*(interv$n[interv$n>min(interv$n[interv$var1>0])]min(interv$n[interv$var1>0]))/(max(interv$n)min(interv$n[interv$var1>0])+1)res.1<-arima(yota_ts, xreg=data.frame(interv$cb,interv$var2),order = c(1,0,1), method="ML")res.1sum(res.1$residuals*res.1$residuals)/length(res.1$residuals)cwp(res.1)fin<rbind(fin,c('[c,b]','complex',res.1$arma,sum(res.1$residuals*res.1$residuals)/length(res.1$residuals),res.1$coef[4],res.1$coef[5]))##################################################################1.
Simple ARIMA#Case CC ready for (1,0,0)(1,0,0)a1<-1a2<-0.6interv<-data.frame(row(as.matrix(data[,4])),data[,3:4])names(interv)<-c('n','var1','var2')interv$cb<-0interv$cb[interv$n>min(interv$n[interv$var1>0])]<-1a1*(interv$n[interv$n>min(interv$n[interv$var1>0])]min(interv$n[interv$var1>0]))/(max(interv$n)min(interv$n[interv$var1>0])+1)interv$bc<-0interv$bc[interv$n>min(interv$n[interv$var2>0])]<-1a2*(interv$n[interv$n>min(interv$n[interv$var2>0])]min(interv$n[interv$var2>0]))/(max(interv$n)min(interv$n[interv$var2>0])+1)res.1<-arima(yota_ts, xreg=data.frame(interv$cb,interv$bc),order = c(1,0,0), method="ML")#res.1#sum(res.1$residuals*res.1$residuals)/length(res.1$residuals)fin<rbind(fin,c('[c,c]','simple',res.1$arma,sum(res.1$residuals*res.1$residuals)/length(res.1$residuals),res.1$coef[3],res.1$coef[4]))##2.
Auto-arima143#CASE (c,c)interv<-data.frame(row(as.matrix(data[,4])),data[,3:4])names(interv)<-c('n','var1','var2')c1.zmin<-1000000000a1<-seq(0.1, 1, 0.1)a2<-seq(0.1, 1, 0.1)for(i in a1) {for(j in a2) {interv$cb<-interv$var1interv$cb[interv$n>min(interv$n[interv$var1>0])]<-1i*(interv$n[interv$n>min(interv$n[interv$var1>0])]min(interv$n[interv$var1>0]))/(max(interv$n)min(interv$n[interv$var1>0])+1)interv$bc<-interv$var2interv$bc[interv$n>min(interv$n[interv$var2>0])]<-1j*(interv$n[interv$n>min(interv$n[interv$var2>0])]min(interv$n[interv$var2>0]))/(max(interv$n)min(interv$n[interv$var2>0])+1)res.1<-auto.arima(yota_ts,xreg=data.frame(interv$cb,interv$bc),max.p = 4, max.q = 4,max.P = 0, max.Q = 0)c1.z<sum(res.1$residuals*res.1$residuals)/length(res.1$residuals)if(c1.z<c1.zmin){c1.zmin<-c1.zc1.a_out<-ic2.a_out<-jmodelmin<-res.1$armamodel<-res.1interv_cc<-data.frame(interv$cb,interv$bc)}}}#zmin#c1.a_out#c2.a_out#modelminwrite.table(interv_cc,file="interv_cc.csv",sep=";",row.names=F,dec=",")fin<rbind(fin,c('[c,c]','auto',modelmin,zmin,model$coef[2],model$coef[3]))##3.
Complex-arima#Case CC ready for (1,0,0)(1,0,0)a1<-1a2<-0.6interv<-data.frame(row(as.matrix(data[,4])),data[,3:4])names(interv)<-c('n','var1','var2')interv$cb<-0interv$cb[interv$n>min(interv$n[interv$var1>0])]<-1a1*(interv$n[interv$n>min(interv$n[interv$var1>0])]144min(interv$n[interv$var1>0]))/(max(interv$n)min(interv$n[interv$var1>0])+1)interv$bc<-0interv$bc[interv$n>min(interv$n[interv$var2>0])]<-1a2*(interv$n[interv$n>min(interv$n[interv$var2>0])]min(interv$n[interv$var2>0]))/(max(interv$n)min(interv$n[interv$var2>0])+1)res.1<-arima(yota_ts, xreg=data.frame(interv$cb,interv$bc),order = c(2,0,2), method="ML")#res.1#sum(res.1$residuals*res.1$residuals)/length(res.1$residuals)cwp(res.1)fin<rbind(fin,c('[c,c]','complex',res.1$arma,sum(res.1$residuals*res.1$residuals)/length(res.1$residuals),res.1$coef[6],res.1$coef[7]))#### FINISHwrite.table(fin,file="fin_matrix_arima_20130318_with_signs.csv", sep=";",dec=",")145Приложение 3.
Реализация алгоритма для оценки одноговнешнего события в программной среде MATLAB%Final edition with loops and parameters% Solve an Input-Output Fitting problem with a Neural Network% Script generated by NFTOOL% Created Fri Oct 05 12:07:35 GMT+04:00 2012%% This script assumes these variables are defined:%%x - input data.%t - target data.neurons = 5; % 6 or 8%n_inputs = 6; % or 12;%intervention_type = 'bb'; %'bc' or 'cb' or 'cc'best_model_counter=0;n_inp = [12]; %6,12interv_types = [1];for n_inputs=n_inpfor intervention_type=interv_typestrue_min_mse=10000;false_min_mse=10000;true_count=0;false_count=0;t_count=0;clear min_perf;clear est_var1;clear est_var2;clear perf_total;clear cars_inputs;clear interventions;n_cols=n_inputs+1;n_rows=neurons;for x=1:1000clear perf;switch n_inputscase 6load lenta_inputs_6.txt;lenta_inputs = lenta_inputs_6;load intervention_6.txt;interventions=intervention_6;load lenta_targets_6.txt;lenta_targets=lenta_targets_6;case 12load lenta_inputs_12.txt;lenta_inputs = lenta_inputs_12;load intervention_12.txt;146interventions=intervention_12;load lenta_targets_12.txt;lenta_targets=lenta_targets_12;endscaled_inputs=mapminmax(lenta_inputs);[scaled_targets,PS]=mapminmax(lenta_targets);inputs = [scaled_inputs;interventions];targets = scaled_targets;% Create a Fitting NetworkhiddenLayerSize = neurons;net = fitnet(hiddenLayerSize);net=configure(net,inputs,targets);net.layers{1}.transferFcn = 'logsig';net.layers{2}.transferFcn = 'purelin';% Choose Input and Output Pre/Post-Processing Functions% For a list of all processing functions type: help nnprocess% Setup Division of Data for Training, Validation, Testing% For a list of all data division functions type: help nndividenet.divideFcn = 'dividerand'; % Divide data randomlynet.divideMode = 'sample'; % Divide up every samplenet.divideParam.trainRatio = 70/100;net.divideParam.valRatio = 10/100;net.divideParam.testRatio = 20/100;% For help on training function 'trainlm' type: help trainlm% For a list of all training functions type: help nntrainnet.trainFcn = 'trainlm'; % Levenberg-Marquardt% Choose a Performance Function% For a list of all performance functions type: helpnnperformancenet.performFcn = 'mse'; % Mean squared error% Choose Plot Functions% For a list of all plot functions type: help nnplotnet.plotFcns = {'plotperform','plottrainstate','ploterrhist',...'plotregression', 'plotfit'};% Train the Networke = targets-net(inputs);perf(1) = mse(e);i=2;147f=0;%% Begin WHILE ?while(f~=1)for nc=1:n_cols,for nr=1:n_rows,if nc~=n_cols && nr==n_rowsnet.IW{1,1}(nr,nc)=0;elseif nr~=n_rows && nc==n_colsnet.IW{1,1}(nr,nc)=0;endendendnet.trainParam.epochs=1;[net,tr] = train(net,inputs,targets);e = targets-net(inputs);perf(i) = mse(e);% Test the Networkoutputs = net(inputs);errors = gsubtract(targets,outputs);performance = perform(net,targets,outputs);% Recalculate Training, Validation and Test PerformancetrainTargets = targets .* tr.trainMask{1};valTargets = targets .* tr.valMask{1};testTargets = targets .* tr.testMask{1};trainPerformance = perform(net,trainTargets,outputs);valPerformance = perform(net,valTargets,outputs);testPerformance = perform(net,testTargets,outputs);if abs(perf(i)-perf(i-1))<10^(-6)f=1;else f=0; endi=i+1;t_count=t_count+1;end%perfa5=net.LW{2,1}(1,n_rows);w5=net.iw{1,1}(n_rows,n_cols);coef = (max(lenta_targets)-min(lenta_targets))/(1-(-1));iw=net.iw{1,1};lw=net.lw{2,1};b1=net.b{1,1};b2=net.b{2,1};b5=net.b{1,1}(n_rows,1);148min_perf(x)=perf(i-1);est_var2(x)=a5 * coef * (1 / (1 + exp(-1*(w5+b5)))-1 / (1 +exp(w5-b5)));if est_var2(x)< 0if min_perf(x)<=true_min_msetrue_min_mse=min_perf(x);true_est_var2=est_var2(x);true_outputs=mapminmax('reverse',outputs,PS);really_mse=sqrt(mse(exp(mapminmax('reverse',outputs,PS))exp(lenta_targets)));true_min_mse_fin=mse(targets-outputs);best_a5=a5;best_w5=w5;best_iw=iw;best_lw=lw;best_b1=b1;best_b2=b2;best_b5=b5;if true_min_mse~=true_min_mse_finbreak;endendtrue_count=true_count+1;elseif min_perf(x)<=false_min_msefalse_min_mse=min_perf(x);false_est_var2=est_var2(x);endfalse_count=false_count+1;endendbest_model_counter=best_model_counter+1;c_number(best_model_counter)=best_model_counterc_neurons(best_model_counter)=neuronsc_n_inputs(best_model_counter)=n_inputsc_intervention_type(best_model_counter)=intervention_typec_really_mse(best_model_counter)=really_msec_true_min_mse(best_model_counter)=true_min_msec_true_min_mse_fin(best_model_counter)=true_min_mse_finc_true_est_var2(best_model_counter)=true_est_var2c_best_a5(best_model_counter)=best_a5c_best_w5(best_model_counter)=best_w5c_best_b5(best_model_counter)=best_b5c_false_min(best_model_counter)=false_min_msec_false_est_var2(best_model_counter)=false_est_var2c_true_count(best_model_counter)=true_count149c_false_count(best_model_counter)=false_countc_t_count(best_model_counter)=t_countendendoutput_matrix=[c_number;c_neurons;c_n_inputs;c_intervention_type;c_really_mse;c_true_min_mse;c_true_min_mse_fin;c_true_est_var2;c_best_a5;c_best_w5;c_best_b5;c_false_min;c_false_est_var2;c_true_count;c_false_count;c_t_count]csvwrite('output_5neuro_12_1000.csv',output_matrix);150Приложение 4.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
