Статистический анализ и моделирование процессов авторегресии и скользящего среднего

Результат представлен на следующем рисунке:

 

Рисунок 6 – Моделирование случайного процесса по модели АРСС(1,3)

Аналогично получим выборки  для моделей СС(0) и АР(3).

6. Оценка моментных функций смоделированного процесса

Найдем оценки моментных функций случайных последовательностей, полученных в п. 5.

 

Таблица 4 – Сравнение моментных функций

Параметры процесса	Исходный процесс	АР(3)		СС(0)		АРСС(1,3)
		Теория	Выборка	Теория	Выборка	Теория	Выборка
Минимум	- 12.089		-20.9079		-21.6642		-19.1330
Максимум	34.1600		43.8839		40.6897		28.3196
Среднее	7.8975	7.8975	7.8326	7.8975	7.9962	7.8975	8.2795
Дисперсия	65.0763	65.0763	63.3375	65.0763	65.5344	65.0763	61.89
СКО	8.0670	8.0670	7.9585	8.0670	8.0953	8.0670	7.8670
Нормированная корреляционная функция
r(0)	1	1	1	1	1	1	1
r(1)	0.9830	0.9830	0.9820	0	-0.0028	0.9830	0.9824
r(2)	0.9481	0.9481	0.9449	0	0.0362	0.9481	0.9473
r(3)	0.9123	0.9123	0.9060	0	0.0228	0.9123	0.9117
r(4)	0.8797	0.8822	0.8718	0	0.0054	0.8778	0.8796
r(5)	0.8489	0.8569	0.8416	0	0.0001	0.8447	0.8493
r(6)	0.8168	0.8336	0.8126	0	0.0267	0.8128	0.8206
r(7)	0.7865	0.8105	0.7838	0	-0.0158	0.7821	0.7937
r(8)	0.7576	0.7875	0.7555	0	-0.0030	0.7525	0.7678
r(9)	0.7286	0.7647	0.7278	0	-0.0073	0.7241	0.7428
r(10)	0.7010	0.7426	0.7006	0	-0.0064	0.6968	0.7187
СКО		0.0048477	0.0002	7.073979	6.9680	0.000078	0.0007

 

Наряду с таблицей представим и график для сравнения результатов  моделирования наилучшей модели с теоретическими предположениями  и исходными данными:

Рисунок 7 – Сравнение корреляционных функций

Заключение

 

В ходе  данного исследования были построены модели АР, СС и смешанные  модели АРСС. Были найдены теоретические параметры моделей. На основе их сравнения выбрана лучшая модель, которая и использовалась для реализации случайного процесса. Установлено, что такая модель довольно точно имитирует исходную случайную последовательность.

Модели авторегрессии  и случайного среднего имеют большую  практическую ценность. На сегодняшний день актуальна задача, когда необходимо смоделировать случайный процессы по уже имеющимся данным в виде выборки. Модели, которые исследовались в данной работе, позволяют успешно ее решать.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Список использованной литературы

 

1. Тараскин, А.Ф. Статистический анализ временных рядов авторегрессии и скользящего среднего: учебное пособие [Текст] // Самара: СГАУ, 1998. – 56с.
2. Тараскин, А.Ф. Статистическое моделирование и метод Монте–Карло: учебное пособие [Текст] // Самара: СГАУ, 1997. – 62с.
3. Храмов, А.Г. Анализ и моделирование процессов АРСС: интернет-ресурс к курсовой работе [Электронный ресурс] // Самара: СГАУ, 2009.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Приложение A Текст программы

//ЗАДАНИЕ 1:

//Оценить моментные  функции случайной последовательности,

//оценить радиус  корреляции, изобразить графически

//оценку нормированной  корреляционной функции

clear()

 

X = fscanfMat('C:\tsp.txt');//Чтение выборки из файла

F = '%16.4f'; // Формат вывода чисел с плавающей точкой

I = '%d'; // Вормат вывода целых чисел

 

function [RX,MX] = Analiz(X)

  MX = mean(X);

  RX_=corr(X,30);

  RX = RX_';

  DX = RX_(1);

  j=30;

  while ((j>0) & (abs(RX(j)/DX)<1/%e))

    j = j-1;

  end;

  rad = j;

  printf("Минимальное значение  : "+F+"\n",min(X));

  printf("Максимальное  значение : "+F+"\n",max(X));

  printf("Выборочное среднее    : "+F+"\n",MX);

  printf("Выборочная дисперсия  : "+F+"\n",DX);

  printf("Выборочное СКО        : "+F+"\n",sqrt(DX));

  printf("Радиус корреляции     :               "+I+"\n",rad);

  printf("Выборочная нормированная  корреляционная функция :\n");

  printf(F+"\n",RX/DX);

  printf("           ...\n           ...\n           ...\n\n")

endfunction;

 

function [RX,m] = first(X)

  [RX,m] = Analiz(X);

  plot2d(0:10,[RX(1:11)/RX(1) zeros(11,1)],axesflag = 5,style = [5 1],leg = "@Normalized correlation function");

endfunction

 

//ЗАДАНИЕ 2:

//Построить модели  АР, СС, АРСС до третьего порядка  включительно

 

//Нахождение альфа  в модели СС

function [alphas,correct] = MA(rX,N)

  RX = rX;

  num = N;

  function [y] = ma(x)

    for i = 0:num,

      y(i+1) = - RX(i+1),

      for j = 0:num-i

        y(i+1) = y(i+1)+x(j+1)*x(j+1+i);

      end;

     end;

  endfunction;

  [alphas,values,info] = fsolve(zeros(1,num+1),ma);

  correct = %T

  for i = 1:N+1

    if (abs(values(i))>0.000001 | info == 4) then

      correct = %F;

      break;

    end;

  end;

endfunction;

 

//Устойчивость модели

function [y] = stable(betas)

  y = %T;

  len = length(betas)

  if abs(betas(len))>1 then

    y = %F;

  elseif (len == 2) then

    if (abs(betas(1))>1-betas(2)) then

      y = %F;

    end;

  elseif (len == 3) then

    if ((abs(betas(1)+betas(3)) > 1-betas(2)) | (abs(betas(2)+betas(1)*betas(3))>abs(1-betas(3)^2))) then

      y = %F;

    end;

  end;

endfunction

 

//Нахождение альфа  и бета в модели АР

function [alpha,betas,is_correct,is_stable] = AR(rX,M)

  RX = rX;

  num = M;

  function [y] = ar(x)

    for i = 0:num,

      y(i+1) = -RX(i+1);

      for j = 1:num,

        y(i+1) = y(i+1)+x(j) * RX(abs(j-i)+1);

      end;

    end;

    y(1) = y(1)+x(num+1)^2;

  endfunction;

  [coef,values,info] = fsolve(zeros(1,num+1),ar);

  is_correct = %T;

  for i = 1:num+1

    if (abs(values(i))>0.000001 | info == 4) then

      is_correct = %F;

      break;

    end;

  end;

  betas = coef(1:M);

  is_stable = %T;

  if (~stable(betas)) then

    is_stable = %F;

  end;

  alpha = coef(M+1);

endfunction

 

//Смешанная корреляционная функция

function R = Mcorr(k, alphas, betas)

  R = alphas(k+1);

  len = min(k, length(betas));

  for j = 1 : len,

    R = R + betas(j) * Mcorr(k - j, alphas);

  end;

endfunction;

 

// Нахождение коэффициентов  бета в общей модели АРСС

function [betas,is_stable] = arma_b(RX, M, N)

  R = zeros(M, 1);

  R_b = zeros(M, M);

  for i = 1:M,

    R(i) = RX(N+1+i);

    for j = 1:M,

      R_b(i, j) = RX(abs(N-j+i)+1);

    end;

  end;

  betas = linsolve(R_b, -R);

  is_stable = %T;

  if (~stable(betas)) then

    is_stable = %F;

  end;

endfunction;

 

//Нахождение коэффициентов  альфа в общей модели АРСС

function [alph,correct] = arma_a(RX, m, n, betas)

  M = m, N = n;

  function r_s = system(alph)

    for k = 0 : N,

      r_s(k+1) = -RX(k+1);

      for i = k : N,

        r_s(k+1) = r_s(k+1) + alph(i+1) * Mcorr(i - k, alph, betas);

      end;

      for j = 1 : M,

        r_s(k+1) = r_s(k+1) + betas(j) * RX(abs(k - j) + 1);

      end;

    end;

  endfunction;

  [alph,values,info] = fsolve([1 : (N+1)], system);

  correct = %T;

  for i = 1:N+1

    if (abs(values(i))>0.000001 | info == 4) then

      correct =%F;

      break;

    end;

  end;

endfunction;

 

//Общая функция  нахождения параметров модели  АРСС

function [alphas,betas,correct,is_stable] = ARMA(RX,M,N)

  if (M == 0) then

    is_stable = %T;

    [alphas,correct] = MA(RX,N);

    betas = [];

  elseif (N == 0) then

    [alphas,betas,correct,is_stable] = AR(RX,M)

  else

    [betas,is_stable] = arma_b(RX,M,N);

    [alphas,correct] = arma_a(RX,M,N,betas);

  end;

endfunction

 

function second(RX)

  for i = 0:3

    for j = 0:3

      [a,b,c,st] = ARMA(RX,i,j),

      printf("Модель АРСС("+I+","+I+") :",i,j);

      if c then

        cor = " существует, ";

        if st then

          cor = cor + "устойчива :\n";

          printf(cor);

          printf(" альфа : ");

          for k = 1:length(a)

            printf(F+" ",a(k));

          end;

          printf("\n");

          if (length(b)) then

            printf(" бета  : ");

            for k = 1:length(b)

              printf(F+" ",b(k));

            end;

            printf("\n");

          end;

        else

          cor = cor + " но не устойчива :\n";

          printf(cor);

        end;

      else

        printf(" не существует\n");

      end;

    end;

  end;

  printf("\n");

endfunction

 

 

//ЗАДАНИЕ 3

//Рассчитать теоретические  нормированные корреляционные функции

//для каждой из  построенных моделей. На основе  сравнения теоретических и выборочных  функций

//выбрать наиболее  адекватную модель из каждого  класса.

//Построить графики  корреляционных функций для трех  наилучших моделей.

 

//Теоретическая  корреляционная функция

function [r] = T_n_corr(RX,alphas,betas,num);

  function R = Tcorr(RX,alphas,betas, k)

    nm = length(alphas)+length(betas)-1;

    k = abs(k);

    if (k > nm) then

      R = 0;

      M = length(betas);

      for j = 1 : M,

        R = R + betas(j) * Tcorr(RX,alphas,betas, k - j);

      end;

    else

      R = RX(k+1);

    end;

  endfunction;

  for i = 0:num-1

    R(i+1) = Tcorr(RX,alphas,betas,i);

  end;

  r = R/R(1);

endfunction

 

function [e] = T_error(rX,r);//Ошибка модели

  e = 0;

  for i = 1:11,

    e = e + (rX(i)-r(i))^2;

  end;

endfunction

 

function [E] = T_errors(RX)//Ошибки моделей

  for i = 1:4

    for j = 1:4

      [a,b,c,s] = ARMA(RX,i-1,j-1);

      if (c & s) then

        r = T_n_corr(RX,a,b,15);

        E(i,j) = T_error(RX/RX(1),r);

      else

        E(i,j) = %inf;

      end;

    end;

  end;

endfunction

 

function [ar_,ma_,arma] = best_models(errors);

  [min_,k] = min(errors(1:4,1));

  ar_ = k-1;

  [min_,k] = min(errors(1,1:4));

  ma_ = k-1;

  [min_,k] = min(errors(2:4,2:4));

  arma = k';

endfunction

 

function [err,ar_,ma_,arma] = third(RX,num)

  err = T_errors(RX);

  [ar_,ma_,arma] = best_models(err);

  printf("Лучшие модели :\n");

  printf("АР("+I+")\n",ar_);

  [a,b,c,st] = ARMA(RX,ar_,0);

  AR_corr = T_n_corr(RX,a,b,num);

  printf("Нормированная крреляция:\n");

  for i = 1:11

    printf(F+"\n",AR_corr(i));

    printf("\n");

  end;

  printf("CC("+I+")\n",ma_);

  [a,b,c,st] = ARMA(RX,0,ma_);

  MA_corr = T_n_corr(RX,a,b,num);

  printf("Нормированная крреляция:\n");

  for i = 1:11

    printf(F+"\n",MA_corr(i));

    printf("\n");

  end;

  printf("АРМА("+I+","+I+")\n\n",arma(1),arma(2));

  [a,b,c,st] = ARMA(RX,arma(1),arma(2));

  ARMA_corr = T_n_corr(RX,a,b,num);

  printf("Нормированная крреляция:\n");

  for i = 1:11

    printf(F+"\n",ARMA_corr(i));

    printf("\n");

  end;

  scf(1);

  plot2d(0:num-1,AR_corr, axesflag = 5,style = 13);

  plot2d(0:num-1,MA_corr, axesflag = 5,style = 6);

  plot2d(0:num-1,ARMA_corr, axesflag = 5,style = 2);

  plot2d(0:num-1,RX(1:num)/RX(1), axesflag = 5,style = 5);

endfunction

 

//ЗАДАНИЕ 4:

//Построить и  изобразить графически параметрическую  оценку 

//спектральной плотности  для трёх наилучших моделей.

function S = ARMA_SPM(alphas,betas,w)

  S1 = 0;S2 = 1;

  for i = 0:length(alphas)-1

    S1 = S1 + alphas(i+1)*exp(%i*w*i);

  end;

  for i = 1:length(betas)

    S2 = S2 - betas(i)*exp(%i*w*i);

  end;

  S = abs(S1/S2)^2;

endfunction

 

function S = initial_SPM(omega, RX)

  S = abs(RX(1));

  for k = 1 : length(RX)-1,

    S = S + 2 * RX(k + 1) * cos(omega * k);

  end;

  if (S<0) then

    S = 0.01;

  end;

endfunction;

 

function SPM(RX,ar_best,ma_best,arma_best);

  w = [0:0.02:%pi];

  len = length(w);

  [a,b,c,st] = ARMA(RX,ar_best,0);

  for i = 1:len

    s_ar(i) = ARMA_SPM(a,b,w(i));

  end;

  s_ar = s_ar/RX(1);

  [a,b,c,st] = ARMA(RX,0,ma_best);

  for i = 1:len

    s_ma(i) = ARMA_SPM(a,b,w(i));

  end;

  s_ma = s_ma/RX(1);

  [a,b,c,st] = ARMA(RX,arma_best(1),arma_best(2));

  for i = 1:len

    s_arma(i) = ARMA_SPM(a,b,w(i));

  end;

  s_arma = s_arma/RX(1);

  for i = 1:len

    s_source(i) = initial_SPM(w(i),RX(1:10)/RX(1));

  end;

  scf(2);

  plot2d(w,s_source, axesflag = 5,style = 5);

  plot2d(w,s_ar, axesflag = 5,style = 13);

  scf(3);

  plot2d(w,s_source, axesflag = 5,style = 5);

  plot2d(w,s_ma, axesflag = 5,style = 3);

  scf(4);

  plot2d(w,s_arma, axesflag = 5,style = 2);

  plot2d(w,s_source, axesflag = 5,style = 5);

endfunction

 

//ЗАДАНИЕ 5:

//Смоделировать  процесс с использованием лучшей  модели

//Результат изобразить  графически

function X = Gauss(n,a,D)     //функция генерирует нормальную  СВ

  X_norm=zeros(n,1),          //по заданным параметрам 

  sigma=sqrt(D),

  for i=1:n,

    sum=0,

    for j=1:12,

      sum =sum+rand(),

    end;

    X_norm(i) = sigma*(sum-6)+a,

  end;

  X=X_norm;

endfunction;

 

function [X] = modeling_ARMA(alphas,betas,MX,num);

  X_ = zeros(num+1000,1);

  ksi = grand(num + 1000, 1, 'nor', 0, 1);

  N = length(alphas) - 1;

  M = length(betas);

  for i = 0:num+999,

    for j =0:N

      if (i-j>=0) then

        X_(i+1) = X_(i+1) + alphas(j+1)*ksi(i-j+1);

      end;

    end;

    for j =1:M

      if (i-j>=0) then

        X_(i+1) = X_(i+1) + betas(j)*X_(i-j+1);

      end;

    end;

  end;

  X = X_(1001:num+1000)+MX;

endfunction;

 

function [ar_model,ma_model,arma_model] = fifth(source,ar_,ma_,arma)

  w = 0:120;

  RX' = corr(source,30);

  MX = mean(source);

  dx = sqrt(cmoment(source,2));

  mx = ones(121,1)*MX;

  dx1 = mx - dx;

  dx2 = mx + dx;

  [a,b,c,st] = ARMA(RX,ar_,0);

  ar_model = modeling_ARMA(a,b,MX,5000)

  [a,b,c,st] = ARMA(RX,0,ma_);

  ma_model = modeling_ARMA(a,b,MX,5000)

  [a,b,c,st] = ARMA(RX,arma(1),arma(2));

  ARMA_corr = T_n_corr(RX,a,b,11);

  arma_model = modeling_ARMA(a,b,MX,5000)

  model_corr = corr(arma_model,11);

  scf(5);

  plot2d(w,source(1:121), axesflag = 5,style = 1);

  plot2d(w,ar_model(1:121), axesflag = 5,style = 13);

  plot2d(w,mx, axesflag = 5,style = 2);

  plot2d(w,dx1, axesflag = 5,style = 5);

  plot2d(w,dx2, axesflag = 5,style = 5);

  scf(6);

  plot2d(w,source(1:121), axesflag = 5,style = 1);