clear ();

// Initial settings

SOURCE_FILE_NAME = 'D: \temp\source. txt'; // File contains input data

FLOAT_FORMAT = '%16.4f'; // Floating point values representation format

INT_FORMAT = '%d'; // Decimal integer values representation format

EPSILON = 1.0E-6; // Precision

MAX_AR_LEVEL = 3;

MAX_MA_LEVEL = 3;

IMITATION_LENGTH = 5000;

// Input data

x = fscanfMat (SOURCE_FILE_NAME); // Sample

n = length (x); // Sample size

// Helpers

function printMat (M, mformat),

[n, m] = size (M);

for i = 1: n,

for j = 1: m,

printf (mformat + "", M (i, j));

end;

printf ("\n");

end;

endfunction;

// 1. Moment functions

// Correlation function esimation

function R = correlation (k, x)

if (k < 0),

k = - k;

end;

R = 0.0;

n = length (x);

meanx = mean (x);

for i = 1: (n-k),

R = R + (x (i) - meanx) * (x (i+k) - meanx);

end;

R = R / (n - k - 1);

endfunction;

// Normalized correlation function

function r = ncorrelation (k, x)

r = correlation (k, x) / correlation (0, x);

endfunction;

// Correlation distance

function T = corrdist (x)

coefficient = 0.01;

T = coefficient * length (x) - 2;

em1 = exp (-1);

while (T >= 0) & (abs (ncorrelation (T, x)) < em1),

T = T - 1;

end;

T = T + 1;

endfunction;

// Drawing normalized correlation function plot

function corrplot ()

p. thickness = 6;

m = 10;

t = [-m: m];

y = zeros (length (t), 1);

for i = 1: length (t),

y (i) = ncorrelation (t (i), x);

end;

plot2d3 (t, y, axesflag=5, style=2);

a = gca ();

p = a. children. children;

p. thickness = 3;

p. mark_mode = "on";

p. mark_size_unit = "point";

p. mark_style = 11;

p. mark_size = 3;

endfunction;

// Main

meanx = mean (x); // Sample mean

svx = variance (x); // Sample variance

m = 10; // Correlation values count

R = zeros (m, 1);

r = zeros (m, 1);

for i = 0: m,

R (i+1) = correlation (i, x);

r (i+1) = ncorrelation (i, x);

end;

Tcorr = corrdist (x);

printf ("Sample mean: " + FLOAT_FORMAT + "\n", meanx);

printf ("Sample variance: " + FLOAT_FORMAT + "\n", svx);

printf ("Correlation function estimation: \n");

printMat (R, FLOAT_FORMAT);

printf ("Normalized correlation function estimation: \n");

printMat (r, FLOAT_FORMAT);

printf ("Correlation distance: " + INT_FORMAT + "\n", Tcorr);

scf (1);

corrplot ();

// 2. Models building

// Autoregression coefficients search

function betas = ar (x, arLevel, maLevel)

R = zeros (2 * arLevel, 1);

for i = (maLevel - arLevel + 1): (maLevel + arLevel),

R (i - maLevel + arLevel) = correlation (i, x);

end;

Rmm = zeros (arLevel, arLevel);

_n = maLevel;

for i = 1: arLevel,

_m = _n;

for j = 1: arLevel,

Rmm (i, j) = R (_m - maLevel + arLevel);

_m = _m - 1;

end;

_n = _n + 1;

end;

Rm = - R (arLevel + 1: 2 * arLevel);

betas = linsolve (Rmm, Rm);

endfunction;

// Mutual Correlation Function

function rrr = mcorrelation (k, alph, betas)

rrr = alph (k+1);

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

for j = 1: len,

rrr = rrr + betas (j) * mcorrelation (k - j, alph);

end;

endfunction;

// Moving average coefficients search

function alphas = ma (x, arLevel, maLevel, betas)

for i = 0: max ([arLevel, maLevel]),

R (i+1) = correlation (i, x);

end;

function zr = syst (alph)

for k = 0: maLevel,

zr (k+1) = - R (k+1);

for i = k: maLevel,

zr (k+1) = zr (k+1) + alph (i+1) * mcorrelation (i - k, alph, betas);

end;

for j = 1: arLevel,

zr (k+1) = zr (k+1) + betas (j) * R (abs (k - j) + 1);

end;

end;

endfunction;

[alphas, values, info] = fsolve ([1: (maLevel+1)], syst);

for i = 1: length (values),

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

alphas (1) = %i;

break;

end;

end;

endfunction;

// Image vector

function s = image (v)

s = %F;

for i = 1: length (v),

if (imag (v (i)) <> 0) then

s = %T;

break;

end;

end;

endfunction;

// Model stability

function s = stable (betas)

p = poly ([pertrans (-betas) 1], "z", "coeff");

z = roots (p);

s = %T;

for i = 1: length (z),

if (abs (z (i)) >= 1) then

s = %F;

break;

end;

end;

endfunction;

// Main

alphas_list = list ();

betas_list = list ();

for i = 0: MAX_AR_LEVEL,

for j = 0: MAX_MA_LEVEL,

betas = ar (x, i, j);

alphas = ma (x, i, j, betas);

alphas_list ($+1) = alphas;

betas_list ($+1) = betas;

printf ("ARMA (" + INT_FORMAT + "," + INT_FORMAT + ") \n", i, j);

if (image (alphas)) then

printf ("Model does not exist. \n");

continue;

end;

if (~stable (betas)) then

printf ("Model exists, but not stable. \n");

continue;

end;

printf ("alpha: \n")

printMat (alphas, FLOAT_FORMAT);

printf ("beta: \n");

printMat (betas, FLOAT_FORMAT);

end;

end;

// 3. Models Analysis

// Theoretical correlation function for ARMA model

function R = theoretical_corr (betas, startR, k)

nm = length (startR) - 1;

k = abs (k);

if (k > nm) then

R = 0;

M = length (betas);

for j = 1: M,

R = R + betas (j) * theoretical_corr (betas, startR, k - j);

end;

else

R = startR (k + 1);

end;

endfunction;

// Normalized theoretical correlation function for ARMA model

function r = norm_theoretical_corr (betas, startR, k)

r = theoretical_corr (betas, startR, k) / theoretical_corr (betas, startR, 0);

endfunction;

// Quadratic Error

function epsilon = quadratic_error (x, y)

epsilon = 0;

m = min (length (x), length (y));

for j = 1: m,

epsilon = epsilon + (x (j) - y (j)) ^2;

end;

endfunction;

function corrplot2 (betas, R, N, M, m_)

p. thickness = 6;

m = m_;

t = [-m: m];

y = zeros (length (t), 1);

for i = 1: length (t),

y (i) = norm_theoretical_corr (betas, R (1: N + M + 1), t (i));

end;

plot2d3 (t, y, axesflag=5, style=2);

a = gca ();

p = a. children. children;

p. thickness = 3;

p. mark_mode = "on";

p. mark_size_unit = "point";

p. mark_style = 11;

p. mark_size = 3;

endfunction;

function eta = imitate (alphas, betas, meanx, count)

defect = 1000;

eta = zeros (count + defect + 1, 1);

ksi = grand (count + defect + 1, 1, 'nor', 0, 1);

N = length (alphas) - 1;

M = length (betas);

for k = 1: count + defect + 1,eta (k) = 0;

for i = 0: N,

if (k - i > 0) then

eta (k) = eta (k) + alphas (i+1) * ksi (k - i);

end;

end;

for j = 1: M,

if (k - j > 0) then

eta (k) = eta (k) + betas (j) * eta (k - j);

end;

end;

end;

eta = eta (defect + 2: count + defect + 1) + meanx;

endfunction;

m = 10; // Analysis Depth

epsilon = zeros (MAX_AR_LEVEL, MAX_MA_LEVEL);

best_ar_eps = %inf;

best_ma_eps = %inf;

best_arma_eps = %inf;

best_ar_alpha = [];

best_ma_alpha = [];

best_arma_alpha = [];

best_ar_beta = [];

best_arma_beta = [];

R = zeros (m+1, 1);

r = zeros (m+1, 1);

for k = 0: m,

R (k+1) = correlation (k, x);

r (k+1) = R (k+1) / R (1);

end;

r_model = zeros (m+1, 1);

for i = 0: MAX_AR_LEVEL,

for j = 0: MAX_MA_LEVEL,

alphas = alphas_list (i * (MAX_MA_LEVEL + 1) + j + 1);

betas = betas_list (i * (MAX_MA_LEVEL + 1) + j + 1);

for k = 0: m,

r_model (k+1) = norm_theoretical_corr (betas, R (1: (i + j + 1)), k);

end;

epsilon (i+1, j+1) = quadratic_error (r, r_model);

if (i == 0) & (epsilon (1, j+1) < best_ma_eps) then

best_ma_eps = epsilon (1, j+1);

best_ma_alpha = alphas;

elseif (j == 0) & (epsilon (i+1, 1) < best_ar_eps) then

best_ar_eps = epsilon (i+1, 1);

best_ar_alpha = alphas;

best_ar_beta = betas;

elseif (epsilon (i+1, j+1) < best_arma_eps) then

best_arma_eps = epsilon (i+1, j+1);

best_arma_alpha = alphas;

best_arma_beta = betas;

end;

end;

end;

printf ("Epsilon: \n");

printMat (epsilon, FLOAT_FORMAT);

printf ("Best models: \nAR (" + INT_FORMAT + "), MA (" + INT_FORMAT + "), ARMA (" + INT_FORMAT + "," + INT_FORMAT + "). \n", length (best_ar_beta), length (best_ma_alpha) - 1, length (best_arma_beta), length (best_arma_alpha) - 1);

scf (2);

corrplot2 (best_ar_beta, R, length (best_ar_beta), 0, m);

scf (3);

corrplot2 ([], R, 0, length (best_ma_alpha) - 1, m);

scf (4);

corrplot2 (best_arma_beta, R, length (best_arma_beta), length (best_arma_alpha) - 1, m);

eta_ar = imitate (best_ar_alpha, best_ar_beta, meanx, IMITATION_LENGTH);

eta_ma = imitate (best_ma_alpha, [], meanx, IMITATION_LENGTH);

eta_arma = imitate (best_arma_alpha, best_arma_beta, meanx, IMITATION_LENGTH);

r_ar = zeros (m+1, 1);

r_ma = zeros (m+1, 1);

r_arma = zeros (m+1, 1);

r_ar_imit = zeros (m+1, 1);

r_ma_imit = zeros (m+1, 1);

r_arma_imit = zeros (m+1, 1);

for k = 0: m,

r_ar (k + 1) = norm_theoretical_corr (best_ar_beta, R (1: length (best_ar_beta) + 1), k);

r_ma (k + 1) = norm_theoretical_corr ([], R (1: length (best_ma_alpha)), k);

r_arma (k + 1) = norm_theoretical_corr (best_arma_beta, R (1: length (best_arma_alpha) + length (best_arma_beta)), k);

r_ar_imit (k + 1) = ncorrelation (k, eta_ar);

r_ma_imit (k + 1) = ncorrelation (k, eta_ma);

r_arma_imit (k + 1) = ncorrelation (k, eta_arma);

end;

scf (5);

plot2d ([0: m], [r r_ar r_ar_imit], style= [5 3 2], axesflag=5, leg="Source@AR (" + string (length (best_ar_beta)) + ") @Imitation");

scf (6);

plot2d ([0: m], [r r_ma r_ma_imit], style= [5 3 2], axesflag=5, leg="Source@MA (" + string (length (best_ma_alpha) - 1) + ") @Imitation");

scf (7);

plot2d ([0: m], [r r_arma r_arma_imit], style= [5 3 2], axesflag=5, leg="Source@ARMA (" + string (length (best_arma_beta)) + "," + string (length (best_arma_alpha) - 1) + ") @Imitation");

// 4. Power Spectral Density

// Power Spectral Density

function Fi = pow_spec_dens_arma (omega, alphas, betas)

s_up = 0;

for k = 0: length (alphas) - 1,s_up = s_up + alphas (k+1) * exp (%i * k * omega);

end;

s_down = 1;

for k = 1: length (betas),

s_down = s_down - betas (k) * exp (%i * k * omega);

end;

Fi = abs (s_up / s_down) ^2;

endfunction;

function Fi = pow_spec_dens (omega, R)

Fi = R (1);

for k = 1: length (R) - 1,Fi = Fi + 2 * R (k + 1) * cos (omega * k);

end;

endfunction;

// Power Spectral Density Plot

function densplot (alphas, betas, R, R_imit)

omega = [0: 0.01: %pi];

len = length (omega);

dens = zeros (len, 1);

dens_source = zeros (len, 1);

dens_imit = zeros (len, 1);

for j = 1: len,

dens (j) = pow_spec_dens_arma (omega (j), alphas, betas) / R (1);

dens_source (j) = pow_spec_dens (omega (j), R) / R (1);

dens_imit (j) = pow_spec_dens (omega (j), R_imit) / R_imit (1);

end;

str = "";

if (length (alphas) == 1) then

str = "AR (" + string (length (betas)) + ")";

elseif (length (betas) == 0) then

str = "MA (" + string (length (alphas) - 1) + ")";

else

str = "ARMA (" + string (length (betas)) + "," + string (length (alphas) - 1) + ")";

end;

plot2d (omega, [dens_source dens dens_imit], style= [5 3 2], axesflag=5, leg="Source@" + str + "@Imitation");

endfunction;

perc = 0.01;

R_50 = zeros (perc * length (x), 1);

R_ar_imit_50 = zeros (perc * length (x), 1);

R_ma_imit_50 = zeros (perc * length (x), 1);

R_arma_imit_50 = zeros (perc * length (x), 1);

for k = 0: length (R_50) - 1,R_50 (k+1) = correlation (k, x);

R_ar_imit_50 (k+1) = correlation (k, eta_ar);

R_ma_imit_50 (k+1) = correlation (k, eta_ma);

R_arma_imit_50 (k+1) = correlation (k, eta_arma);

end;

scf (8);

densplot (best_ar_alpha, best_ar_beta, R_50, R_ar_imit_50);

scf (9);

densplot (best_ma_alpha, [], R_50, R_ma_imit_50);

scf (10);

densplot (best_arma_alpha, best_arma_beta, R_50, R_arma_imit_50);

// 5. Imitation

function imitation_plot (x, imitation, meanx, sv)

t = [1: 100];

q = [x (t) imitation (t) (zeros (length (t), 1) + meanx) (zeros (length (t), 1) + meanx + sqrt (sv)) (zeros (length (t), 1) + meanx - sqrt (sv))];

plot2d (t, q, style= [5 2 4 3 3], axesflag=5, leg="Source@Imitation@Mean@Standard deviation");

endfunction;

best_alphas = [];

best_betas = [];

if (best_ar_eps < best_ma_eps) then

if (best_ar_eps < best_arma_eps) then

best_alphas = best_ar_alpha;

best_betas = best_ar_beta;

else

best_alphas = best_arma_alpha;

best_betas = best_arma_beta;

end;

else

if (best_ma_eps < best_arma_eps) then

best_alphas = best_ma_alpha;

best_betas = best_ma_beta;

else

best_alphas = best_arma_alpha;

best_betas = best_arma_beta;

end;

end;

printf ("Best model: ARMA (" + INT_FORMAT + "," + INT_FORMAT + ") \n", length (best_betas), length (best_alphas) - 1);

for i = 1: length (r_model),

r_model (i) = norm_theoretical_corr (best_betas, R (1: (length (best_betas) + length (best_alphas))), i-1);

end;

imitation = imitate (best_alphas, best_betas, meanx, IMITATION_LENGTH);

scf (11);

imitation_plot (x, imitation, meanx, svx);

// 6. Imitation Analysis

function total_sample_analysis (x)

m = 10;

R = zeros (m + 1, 1);

R_src = zeros (m + 1, 1);

for k = 0: m,

R (k + 1) = correlation (k, x);

end;

printf ("Minimum: " + FLOAT_FORMAT + "\n", min (x));

printf ("Maximum: " + FLOAT_FORMAT + "\n", max (x));

printf ("Mean: " + FLOAT_FORMAT + "\n", mean (x));

printf ("Variance: " + FLOAT_FORMAT + "\n", R (1));

printf ("Standard devation: " + FLOAT_FORMAT + "\n", sqrt (R (1)));

printf ("Normalized correlation function: \n");

printMat (R / R (1), FLOAT_FORMAT);

printf ("Epsilon: " + FLOAT_FORMAT + "\n", quadratic_error (R / R (1), r));

endfunction;

function total_model_analysis (alphas, betas, R)

R_model = zeros (length (R), 1);

for k = 0: m,

R_model (k + 1) = theoretical_corr (betas, R (1: length (alphas) + length (betas)), k);

end;

printf ("Mean: " + FLOAT_FORMAT + "\n", meanx);

printf ("Variance: " + FLOAT_FORMAT + "\n", R_model (1));

printf ("Standard devation: " + FLOAT_FORMAT + "\n", sqrt (R_model (1)));

printf ("Normalized correlation function: \n");

printMat (R_model / R_model (1), FLOAT_FORMAT);

printf ("Epsilon: " + FLOAT_FORMAT + "\n", quadratic_error (R_model / R_model (1), r));

endfunction;

mean_imit = mean (imitation);

var_imit = variance (imitation);

m = 10;

r_imit = zeros (m, 1);

for k = 0: m,

r_imit (k+1) = ncorrelation (k, imitation);

end;

printf ("Imitation mean: " + FLOAT_FORMAT + "\n", mean_imit);

printf ("Imitation variance: " + FLOAT_FORMAT + "\n", var_imit);

printf ("Imitation normalized correlation function: \n");

printMat (r_imit, FLOAT_FORMAT);

printf ("ARMA model theoretical correlation function: \n");

printMat (r_model, FLOAT_FORMAT);

printf ("\nTotal table\n");

printf ("\nSource\n");

total_sample_analysis (x);

printf ("\nAR (" + INT_FORMAT + ") \n", length (best_ar_beta));

printf ("Theory\n");

total_model_analysis (best_ar_alpha, best_ar_beta, R);

printf ("Sample\n");

total_sample_analysis (eta_ar);

printf ("\nMA (" + INT_FORMAT + ") \n", length (best_ma_alpha) - 1);

printf ("Theory\n");

total_model_analysis (best_ma_alpha, [], R);

printf ("Sample\n");

total_sample_analysis (eta_ma);

printf ("\nARMA (" + INT_FORMAT + "," + INT_FORMAT + ") \n", length (best_arma_beta), length (best_arma_alpha) - 1);

printf ("Theory\n");

total_model_analysis (best_arma_alpha, best_arma_beta, R);

printf ("Sample\n");

total_sample_analysis (eta_arma);

Размещено на Allbest.ru