Документация

Модель и аналитический обзор

Цена электричества сопоставлена с соответствующим спросом [1]. Изменения в численности населения и техническом прогрессе предлагают, чтобы спотовые цены электричества имели долгосрочный тренд. Кроме того, учитывая климат области, спрос на электричество колеблется согласно сезону. Следовательно, спотовые цены электричества колеблются, точно так же предлагая включение сезонного компонента в модели.

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

Характеристики спотовых цен электричества влияют на шаги для образцового создания:

Загрузите и визуализируйте спотовые цены электричества

Загрузите MAT-файл Data_ElectricityPrices, включенный с Econometrics Toolbox™. Данные состоят из 1411 1 расписание DataTable, содержащий ежедневные спотовые цены электричества.

load Data_ElectricityPrices
T = size(DataTable,1); % Sample size

Рабочая область содержит 1411 1 расписание MATLAB® DataTable ежедневных спотовых цен электричества среди других переменных.

Определите, содержат ли данные тренды путем графического вывода спотовых цен в зависимости от времени.

figure
plot(DataTable.Date, DataTable.SpotPrice)
xlabel('Date')
ylabel('Spot price ($)')
title('Daily Electricity Prices')
grid on
axis tight

Выставка спотовых цен:

Присутствие экспоненциального тренда трудно определить из графика.

Оцените и удалите экспоненциальный тренд

Оцените присутствие экспоненциального тренда путем вычисления ежемесячной статистики. Если ряд спотовой цены показывает экспоненциальный тренд, то средние значения и стандартные отклонения, вычисленные по постоянным периодам времени, лежат на строке со значительным наклоном.

Оцените ежемесячные средние значения и стандартные отклонения ряда.

statsfun = @(v) [mean(v) std(v)];
MonthlyStats = retime(DataTable,'monthly',statsfun);

Постройте ежемесячные стандартные отклонения по средним значениям. Добавьте строку лучшей подгонки к графику.

figure
scatter(MonthlyStats.SpotPrice(:, 1),...
        MonthlyStats.SpotPrice(:, 2),'filled')
h = lsline;
h.LineWidth = 2.5;
h.Color = 'r';
xlabel('Mean ($)')
ylabel('Standard deviation ($)')
title('Monthly Price Statistics')
grid on

Выполните тест для значения линейной корреляции. Включайте результаты в легенду.

[rho, p] = corr(MonthlyStats.SpotPrice);
legend({'(\mu, \sigma)', ...
    ['Line of best fit (\rho = ',num2str(rho(1, 2),'%0.3g'),...
    ', {\it p}-value: ',num2str(p(1, 2),'%0.3g'),')']},...
    'Location','northwest')

P-значение близко к 0, предполагая, что ежемесячные статистические данные значительно положительно коррелируются. Этот результат представляет свидетельства экспоненциального тренда в ряду спотовой цены.

Удалите экспоненциальный тренд путем применяния логарифмического преобразования к ряду.

DataTable.LogPrice = log(DataTable.SpotPrice);

Стройте цены журнала в зависимости от времени.

figure
plot(DataTable.Date,DataTable.LogPrice)
xlabel('Date')
ylabel('Log price')
title('Daily Log Electricity Prices')
grid on
axis tight

Оцените долгосрочный детерминированный тренд

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

ElapsedYears = years(DataTable.Date - DataTable.Date(1)); % Number of elapsed years
DataTable = addvars(DataTable,ElapsedYears,'Before',1);   % ElapsedYears is first variable in DataTable

lccoeffs = polyfit(DataTable.ElapsedYears,DataTable.LogPrice,1);
DataTable.LinearTrend = polyval(lccoeffs,DataTable.ElapsedYears);

Постройте линейный компонент тренда по логарифмическим спотовым ценам.

hold on
plot(DataTable.Date,DataTable.LinearTrend,'LineWidth',2)
legend({'Log price' 'Linear trend'})

Идентифицируйте сезонные компоненты долгосрочного тренда путем выполнения спектрального анализа данных. Запустите анализ путем выполнения этих действий:

DataTable.LinearDetrendedLogPrice = DataTable.LogPrice - DataTable.LinearTrend;
fftLogPrice = fft(DataTable.LinearDetrendedLogPrice);

powerLogPrice = fftLogPrice.*conj(fftLogPrice)/T; % Power spectrum
freq = (0:T - 1)*(1/T);                           % Frequency vector

Поскольку спектр сигнала с действительным знаком симметричен о средней частоте, удалите последнюю половину наблюдений от данных о частоте и степени.

m = ceil(T/2);
powerLogPrice = powerLogPrice(1:m);
freq = freq(1:m);

Преобразуйте частоты в периоды. Поскольку наблюдения ежедневно измеряются, делят периоды на 365,25, чтобы выразить их в годах.

periods = 1./freq/365.25;

Идентифицируйте два периода, которые имеют самые большие степени.

[maxPowers,posMax] = maxk(powerLogPrice,2);
topPeriods = periods(posMax);

figure;
stem(periods,powerLogPrice,'.')
xlabel('Period (years)')
ylabel('Power')
title('Power Spectrum of Daily Log Electricity Prices')
hold on
plot(periods(posMax),maxPowers,'r^','MarkerFaceColor','r')
grid on
legend({'Power','Maxima'})

Доминирующие сезонные компоненты соответствуют 6-месячным и 12-месячным циклам.

Соответствуйте долгосрочной детерминированной модели тренда

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

где:

Синус и члены в модели косинуса составляют возможное смещение фазы. Таким образом, поскольку фаза смещается $$\theta$$:

$$A\sin (ft+\theta )=(A\mathrm{потому\; что}\theta )\sin ft+(A\sin \theta )\mathrm{потому\; что}ft=B\sin ft+C\mathrm{потому\; что}ft$$,

где $$B=A\mathrm{потому\; что}\theta$$ и $$C=A\sin \theta$$ константы. Поэтому включение синуса и условий косинуса с той же частотой составляет смещение фазы.

Сформируйте матрицу проекта. Поскольку программное обеспечение включает постоянный член в модели по умолчанию, не включайте столбец из единиц в матрице проекта.

designMat = @(t) [t sin(2*pi*t) cos(2*pi*t) sin(4*pi*t) cos(4*pi*t)];
X = designMat(DataTable.ElapsedYears);

Соответствуйте модели к логарифмическим спотовым ценам. Примените пошаговую регрессию, чтобы выбрать важные предикторы.

varNames = {'t' 'sin(2*pi*t)' 'cos(2*pi*t)'...
            'sin(4*pi*t)' 'cos(4*pi*t)' 'LogPrice'};
TrendMdl = stepwiselm(X,DataTable.LogPrice,'Upper','linear','VarNames',varNames);

1. Adding t, FStat = 109.7667, pValue = 8.737506e-25
2. Adding cos(2*pi*t), FStat = 135.2363, pValue = 6.500039e-30
3. Adding cos(4*pi*t), FStat = 98.0171, pValue = 2.19294e-22
4. Adding sin(4*pi*t), FStat = 22.6328, pValue = 2.16294e-06

disp(TrendMdl)

Linear regression model:
    LogPrice ~ 1 + t + cos(2*pi*t) + sin(4*pi*t) + cos(4*pi*t)

Estimated Coefficients:
                   Estimate        SE         tStat       pValue  
                   _________    _________    _______    __________

    (Intercept)       3.8617     0.014114      273.6             0
    t              -0.073594    0.0063478    -11.594    9.5312e-30
    cos(2*pi*t)     -0.11982     0.010116    -11.845    6.4192e-31
    sin(4*pi*t)     0.047563    0.0099977     4.7574    2.1629e-06
    cos(4*pi*t)     0.098425    0.0099653     9.8768    2.7356e-22


Number of observations: 1411, Error degrees of freedom: 1406
Root Mean Squared Error: 0.264
R-squared: 0.221,  Adjusted R-Squared: 0.219
F-statistic vs. constant model: 99.9, p-value = 7.15e-75

stepwiselm понижается $$\sin (2\pi t)$$ назовите из линейной модели, указав, что ежегодный цикл начинается на пике волны.

Постройте подобранную модель.

DataTable.DeterministicTrend = TrendMdl.Fitted;

figure
plot(DataTable.Date,DataTable.LogPrice)
xlabel('Date')
ylabel('Log price')
title('Daily Log Electricity Prices')
grid on
hold on
plot(DataTable.Date,DataTable.DeterministicTrend,'LineWidth',2)
legend({'Log price','Deterministic trend'})
axis tight
hold off

Анализируйте детрендированные данные

Сформируйте детрендированные временные ряды путем удаления предполагаемого долгосрочного детерминированного тренда из логарифмических спотовых цен. Другими словами, извлеките невязки из подходящей пошаговой модели линейной регрессии.

DataTable.DetrendedLogPrice = TrendMdl.Residuals.Raw;

Отображайте детрендированные данные на графике в зависимости от времени.

figure
plot(DataTable.Date,DataTable.DetrendedLogPrice)
xlabel('Date')
ylabel('Residual')
title('Trend Model Residuals')
grid on
axis tight

Детрендированные данные кажутся в центре в нуле, и ряд показывает последовательную автокорреляцию, потому что несколько кластеров последовательных невязок происходят выше и ниже $\mathit{y}=0$. Эти функции предполагают, что авторегрессивная модель подходит для детрендированных данных.

Чтобы определить количество задержек, чтобы включать в авторегрессивную модель, примените методологию Поля-Jenkins. Постройте автокорреляционную функцию (ACF) и частичная автокорреляционная функция (PACF) детрендированных данных в той же фигуре, но на отдельных графиках. Исследуйте первые 50 задержек.

numLags = 50;

figure
subplot(2,1,1)
autocorr(DataTable.DetrendedLogPrice,numLags)
subplot(2,1,2)
parcorr(DataTable.DetrendedLogPrice,numLags)

ACF постепенно затухает. PACF показывает резкое сокращение после первой задержки. Это поведение показательно из AR (1) процесс с формой

$${\xi }_t=\varphi {\xi }_{t-1}+u_t$$,

где $$\varphi$$ задержка 1 авторегрессивный коэффициент и $$u_t\sim N(0,{\sigma }^2)$$ iid последовательность случайных переменных. Поскольку детрендированные данные сосредоточены вокруг нуля, образцовая константа 0.

Создайте модель AR (1) для детрендированных данных без образцовой константы.

DTMdl = arima('Constant',0,'ARLags',1);
disp(DTMdl)

  arima with properties:

     Description: "ARIMA(1,0,0) Model (Gaussian Distribution)"
    Distribution: Name = "Gaussian"
               P: 1
               D: 0
               Q: 0
        Constant: 0
              AR: {NaN} at lag [1]
             SAR: {}
              MA: {}
             SMA: {}
     Seasonality: 0
            Beta: [1×0]
        Variance: NaN

DTMdl является объектом модели arima, представляющим модель AR (1) и служащим шаблоном для оценки. Свойства DTMdl со значением NaN соответствуют допускающим оценку параметрам в модели AR (1).

Соответствуйте модели AR (1) к детрендированным данным.

EstDTMdl = estimate(DTMdl,DataTable.DetrendedLogPrice);

 
    ARIMA(1,0,0) Model (Gaussian Distribution):
 
                 Value      StandardError    TStatistic      PValue   
                ________    _____________    __________    ___________

    Constant           0              0           NaN              NaN
    AR{1}         0.4818       0.024353        19.784       4.0787e-87
    Variance    0.053497      0.0014532        36.812      1.1867e-296

estimate приспосабливает все допускающие оценку параметры в DTMdl к данным, затем возвращает EstDTMdl, объект модели arima, представляющий подходящую модель AR (1).

Выведите невязки и условное отклонение из подходящей модели AR (1).

[DataTable.Residuals,condVar] = infer(EstDTMdl,DataTable.DetrendedLogPrice);

condVar является T-by-1 вектор условных отклонений для каждого момента времени. Поскольку заданная модель ARIMA гомоскедастична, все элементы condVar равны.

Стандартизируйте невязки путем деления каждой невязки на мгновенное стандартное отклонение.

DataTable.StandardizedResiduals = DataTable.Residuals./sqrt(condVar);

Постройте стандартизированные невязки и проверьте, что невязки не автокоррелируются путем графического вывода их ACF к 50 задержкам.

figure
subplot(2,1,1)
plot(DataTable.Date,DataTable.StandardizedResiduals)
xlabel('Date')
ylabel('Residual')
title('Standardized Residuals')
grid on
axis tight
subplot(2,1,2)
autocorr(DataTable.StandardizedResiduals,numLags)

Невязки не показывают автокорреляцию.

Постройте гистограмму стандартизированных невязок и наложите непараметрическую подгонку ядра.

kernelFit = fitdist(DataTable.StandardizedResiduals,'kernel');
sampleRes = linspace(min(DataTable.StandardizedResiduals),...
    max(DataTable.StandardizedResiduals),500).';
kernelPDF = pdf(kernelFit,sampleRes);

figure
histogram(DataTable.StandardizedResiduals,'Normalization','pdf')
xlabel('Residual')
ylabel('Density')
s = skewness(DataTable.StandardizedResiduals);
title(sprintf('\\bf Residual Distribution (Skewness $\\gamma$ = %4f)',s),'Interpreter','latex')
grid on
hold on
plot(sampleRes,kernelPDF,'LineWidth',2)
legend({'Standardized Residuals','Kernel Fit'})
hold off

Создайте график нормального распределения стандартизированных невязок.

figure
probplot(DataTable.StandardizedResiduals)
grid on

Невязки показывают положительную скошенность, потому что они отклоняются от нормальности в верхнем хвосте.

Образцовый скошенный остаточный ряд

Скошенное нормальное распределение эпсилона является семейством почти нормальных распределений с местоположением $\mu$шкала $\sigma$, и дополнительный параметр скошенности $\epsilon$. Модели параметра скошенности любая ненулевая скошенность в данных [2]. Если $$\epsilon =0$$, скошенное нормальное распределение эпсилона уменьшает до нормального распределения.

Econometrics Toolbox™ включает объект пользовательского дистрибутива представление скошенного нормального распределения эпсилона в mlr/examples/econ/helper/+prob/EpsilonSkewNormalDistribution.m, где mlr является значением matlabroot. Чтобы создать ваш собственный объект распределения, выполните эти шаги:

Чтобы гарантировать, что среда распределения вероятностей распознает новое распределение, сбрасывает его список рассылки.

makedist -reset

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

ResDist = fitdist(DataTable.StandardizedResiduals,'EpsilonSkewNormal');
disp(ResDist)

  EpsilonSkewNormalDistribution

  EpsilonSkewNormal distribution
      theta = -0.421946   [-0.546991, -0.296901]
      sigma =  0.972487   [0.902701, 1.04227]
    epsilon = -0.286248   [-0.360485, -0.212011]

Предполагаемый параметр скошенности $$\underset{}{\overset{ˆ}{\epsilon }}$$ приблизительно-0.286, указывая, что невязки положительно скашиваются.

Отобразите предполагаемые стандартные погрешности оценок параметра.

stdErr = sqrt(diag(ResDist.ParameterCovariance));
SETbl = table(ResDist.ParameterNames',stdErr,...
    'VariableNames',{'Parameter', 'StandardError'});

disp('Estimated parameter standard errors:')

Estimated parameter standard errors:

disp(SETbl)

    Parameter    StandardError
    _________    _____________

    'theta'          0.0638   
    'sigma'        0.035606   
    'epsilon'      0.037877   

$$\underset{}{\overset{ˆ}{\epsilon }}$$ имеет небольшую стандартную погрешность (~0.038), указывая, что скошенность является значительной.

Постройте гистограмму невязок и наложите подходящее распределение.

fitVals = pdf(ResDist,sampleRes);

figure
histogram(DataTable.StandardizedResiduals,'Normalization','pdf')
xlabel('Residual')
ylabel('Density')
title('Residual Distribution')
grid on
hold on
plot(sampleRes,fitVals,'LineWidth',2)
legend({'Residuals','Epsilon-Skew-Normal Fit'})
hold off

Подходящее распределение, кажется, хорошая модель для невязок.

Оцените качество подгонки

Оцените качество подгонки скошенного нормального распределения эпсилона путем выполнения этой процедуры:

Вычислите эмпирический cdf стандартизированных невязок и cdf подходящего скошенного нормального распределения эпсилона. Возвратите точки, в которых программное обеспечение оценивает эмпирический cdf.

[resCDF,queryRes] = ecdf(DataTable.StandardizedResiduals);
fitCDF = cdf(ResDist,sampleRes);

Постройте эмпирический и приспособленный cdfs в той же фигуре.

figure
stairs(queryRes,resCDF,'LineWidth', 1.5)
hold on
plot(sampleRes,fitCDF,'LineWidth', 1.5)
xlabel('Residual')
ylabel('Cumulative probability')
title('Empirical and Fitted CDFs (Standardized Residuals)')
grid on
legend({'Empirical CDF','Epsilon-Skew-Normal CDF'},...
        'Location','southeast')

Эмпирические и ссылочные cdfs кажутся подобными.

Проведите тест Кольмогорова-Смирнова для качества подгонки путем выполнения этой процедуры:

rng default                          % For reproducibility
cvp = cvpartition(T,'HoldOut',0.20); % Define data partition on T observations
idxtraining = training(cvp);         % Extract training set indices
idxtest = test(cvp);                 % Extract test set indices

trainingset = DataTable.StandardizedResiduals(idxtraining);
testset = DataTable.StandardizedResiduals(idxtest);    
resFitTrain = fitdist(trainingset,'EpsilonSkewNormal');     % Reference distribution

[~,kspval] = kstest(testset,'CDF',resFitTrain);
disp(['Kolmogorov-Smirnov test p-value: ',num2str(kspval)])

Kolmogorov-Smirnov test p-value: 0.85439

P-значение теста является достаточно большим, чтобы предложить, чтобы нулевая гипотеза, что дистрибутивы являются тем же самым, не была отклонена.

Постройте максимум cdf несоответствие.

fitCDF = cdf(ResDist,queryRes);  % Fitted cdf with respect to empirical cdf query points         
[maxDiscrep,maxPos] = max(abs(resCDF - fitCDF));
resAtMaxDiff = queryRes(maxPos);
plot(resAtMaxDiff*ones(1, 2),...
    [resCDF(maxPos) fitCDF(maxPos)],'k','LineWidth',2,...
    'DisplayName',['Maximum Discrepancy (%): ',num2str(100*maxDiscrep,'%.2f')])

Максимальное несоответствие является маленьким, предлагая, чтобы стандартизированные невязки следовали за заданным скошенным нормальным распределением эпсилона.

Моделируйте будущие спотовые цены электричества

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

Сохраните данные моделирования в расписании под названием SimTbl.

numSteps = 2 * 365;
Date = DataTable.Date(end) + days(1:numSteps).';
SimTbl = timetable(Date);

Моделируйте 1 000 путей Монте-Карло стандартизированных невязок от подходящего скошенного нормального распределения эпсилона ResDist.

numPaths = 1000;
SimTbl.StandardizedResiduals = random(ResDist,[numSteps numPaths]); 

SimTbl.StandardizedResiduals является numSteps-by-numPaths матрица моделируемых стандартизированных остаточных путей. Строки соответствуют периодам в горизонте прогноза, и столбцы соответствуют отдельным путям к симуляции.

Моделируйте 1 000 путей детрендированных цен журнала путем пропущения моделируемых, стандартизированных невязок через предполагаемую модель AR (1) EstDTMdl. Задайте предполагаемые условные отклонения condVar как преддемонстрационные условные отклонения, наблюдаемые детрендированные цены журнала в DataTable как преддемонстрационные ответы и предполагаемые стандартизированные невязки в DataTable как преддемонстрационные воздействия. В этом контексте "предварительная выборка" является периодом перед горизонтом прогноза.

SimTbl.DetrendedLogPrice = filter(EstDTMdl,SimTbl.StandardizedResiduals, ...
    'Y0', DataTable.DetrendedLogPrice,'V0',condVar,...
    'Z0', DataTable.StandardizedResiduals);

SimTbl.DetrendedLogPrice является матрицей детрендированных цен журнала с теми же размерностями как SimTbl.StandardizedResiduals.

Предскажите, что долгосрочный детерминированный тренд путем оценки предполагаемой линейной модели TrendMdl в номере лет протек с начала выборки.

SimTbl.ElapsedYears = years(SimTbl.Date - DataTable.Date(1));
SimTbl.DeterministicTrend = predict(TrendMdl, designMat(SimTbl.ElapsedYears));

SimTbl.DeterministicTrend является долгосрочным детерминированным компонентом тренда логарифмических спотовых цен в горизонте прогноза.

Моделируйте спотовые цены backtransforming моделируемые значения. Таким образом, объедините предсказанный долгосрочный детерминированный тренд и моделируемые логарифмические спотовые цены, затем exponentiate сумма.

SimTbl.LogPrice = SimTbl.DetrendedLogPrice + SimTbl.DeterministicTrend;
SimTbl.SpotPrice = exp(SimTbl.LogPrice);

SimTbl.SpotPrice является матрицей моделируемых путей к спотовой цене с теми же размерностями как SimTbl.StandardizedResiduals.

Анализируйте моделируемые пути

Постройте моделируемый путь представителя. Представительный путь является путем, который имеет окончательное значение, самое близкое к среднему окончательному значению среди всех путей.

Извлеките моделируемые значения в конце горизонта прогноза, затем идентифицируйте путь, который является самым близким к медиане окончательных значений.

finalValues = SimTbl.SpotPrice(end, :);
[~,pathInd] = min(abs(finalValues - median(finalValues)));

Отобразите на графике исторические данные и детерминированный тренд и выбранный моделируемый путь, и предсказал детерминированный тренд.

figure
plot(DataTable.Date,DataTable.SpotPrice)
hold on
plot(DataTable.Date,exp(DataTable.DeterministicTrend),'LineWidth',2)
plot(SimTbl.Date,SimTbl.SpotPrice(:,pathInd),'Color',[1, 0.5, 0])
plot(SimTbl.Date,exp(SimTbl.DeterministicTrend),'k','LineWidth', 2)
xlabel('Date')
ylabel('Spot price ($)')
title('Electricity Spot Prices')
legend({'Historical spot price','Historical trend',...
    'Simulated spot price','Simulated trend'})
grid on
hold off

Получите статистику Монте-Карло из моделируемых путей к спотовой цене путем вычисления, для каждого момента времени в горизонте прогноза, среднем значении и медиане, и 2.5th, 5-й, 25-й, 75-й, 95-й, и 97.5th процентили.

SimTbl.MCMean = mean(SimTbl.SpotPrice,2);
SimTbl.MCQuantiles = quantile(SimTbl.SpotPrice,[0.025 0.05 0.25 0.5 0.75 0.95 0.975],2);

Постройте исторические спотовые цены и детерминированный тренд, все предсказанные пути, статистику Монте-Карло и предсказанный детерминированный тренд.

figure
hHSP = plot(DataTable.Date,DataTable.SpotPrice);
hold on
hHDT = plot(DataTable.Date,exp(DataTable.DeterministicTrend),'LineWidth', 2);
hFSP = plot(SimTbl.Date,SimTbl.SpotPrice,'Color',[0.9,0.9,0.9]);
hFDT = plot(SimTbl.Date,exp(SimTbl.DeterministicTrend),'y','LineWidth', 2);
hFQ = plot(SimTbl.Date,SimTbl.MCQuantiles,'Color',[0.7 0.7 0.7]);
hFM = plot(SimTbl.Date,SimTbl.MCMean,'r--');
xlabel('Date')
ylabel('Spot price ($)')
title('Electricity Spot Prices')
h = [hHSP hHDT hFSP(1) hFDT hFQ(1) hFM];
legend(h,{'Historical spot price','Historical trend','Forecasted paths',...
    'Forecasted trend','Monte Carlo quantiles','Monte Carlo mean'})
grid on

Также, если у вас есть лицензия Financial Toolbox™, можно использовать fanplot, чтобы построить квантили Монте-Карло.

Оцените, обращается ли модель к большим скачкам, показанным в исторических данных путем выполнения этой процедуры:

Вычислите ежемесячные средние значения Монте-Карло и стандартные отклонения моделируемых путей к спотовой цене.

SimSpotPrice = timetable(SimTbl.Date, SimTbl.SpotPrice(:, end),...
    'VariableNames',{'SpotPrice'});
SimMonthlyStats = retime(SimSpotPrice,'monthly',statsfun);

Постройте исторические ежемесячные моменты и их строку лучшей подгонки. Наложите моменты Монте-Карло.

figure
hHMS = scatter(MonthlyStats.SpotPrice(:, 1),...
        MonthlyStats.SpotPrice(:, 2),'filled');
hHL = lsline;
hHL.LineWidth = 2.5;
hHL.Color = hHMS.CData;
hold on
scatter(SimMonthlyStats.SpotPrice(:, 1),...
        SimMonthlyStats.SpotPrice(:, 2),'filled');

Оцените и постройте полный график лучшей подгонки.

mmn = [MonthlyStats.SpotPrice(:, 1); SimMonthlyStats.SpotPrice(:, 1)];
msd = [MonthlyStats.SpotPrice(:, 2); SimMonthlyStats.SpotPrice(:, 2)];
p = polyfit(mmn,msd,1);
overallTrend = polyval(p,mmn);
plot(mmn,overallTrend,'k','LineWidth',2.5)    
xlabel('Monthly mean price ($)')
ylabel('Monthly price standard deviation ($)')
title('Monthly Price Statistics')
legend({'Historical moments','Historical trend',...
    'Simulated moments','Overall trend'},'Location','northwest')
grid on
hold off

Моделируемые ежемесячные моменты широко сопоставимы с историческими данными. Моделируемые спотовые цены имеют тенденцию показывать меньшие скачки, чем наблюдаемые спотовые цены. Чтобы составлять большие скачки, можно смоделировать больший правый хвост путем применения теории экстремума на шаге подбора кривой распределения. Этот подход использует скошенное нормальное распределение эпсилона для невязок, но моделирует верхний хвост при помощи обобщенного распределения Парето. Для получения дополнительной информации смотрите Используя Теорию Экстремума и Связки, чтобы Оценить Риск Рынка.

Сравните скос эпсилона нормальные результаты с предположением нормального распределения

В предыдущем анализе стандартизированные невязки выводят от подходящего скошенного нормального распределения эпсилона. Получите моделируемые спотовые цены backtransforming нормально распределенными, стандартизированными невязками. Можно моделировать стандартизированные невязки непосредственно при помощи функции simulate объектов модели arima. Сохраните компоненты временных рядов в таблице под названием SimNormTbl.

SimNormTbl = timetable(Date);
SimNormTbl.ElapsedYears = SimTbl.ElapsedYears;

SimNormTbl.DeterministicTrend = SimTbl.DeterministicTrend;

SimNormTbl.DetrendedLogPrice = simulate(EstDTMdl,numSteps,'NumPaths',numPaths,...
    'E0',DataTable.Residuals,'V0',condVar,'Y0',DataTable.DetrendedLogPrice);
SimNormTbl.LogPrice = SimNormTbl.DetrendedLogPrice + SimNormTbl.DeterministicTrend;
SimNormTbl.SpotPrice = exp(SimNormTbl.LogPrice);

Получите статистику Монте-Карло из моделируемых путей к спотовой цене путем вычисления, для каждого момента времени в горизонте прогноза, среднем значении и медиане, и 2.5th, 5-й, 25-й, 75-й, 95-й, и 97.5th процентили.

SimNormTbl.MCMean = mean(SimNormTbl.SpotPrice,2);
SimNormTbl.MCQuantiles = quantile(SimNormTbl.SpotPrice,[0.025 0.05 0.25 0.5 0.75 0.95 0.975],2);

Постройте исторические спотовые цены и детерминированный тренд, все предсказанные пути, статистику Монте-Карло и предсказанный детерминированный тренд.

figure
hHSP = plot(DataTable.Date,DataTable.SpotPrice);
hold on
hHDT = plot(DataTable.Date,exp(DataTable.DeterministicTrend),'LineWidth', 2);
hFSP = plot(SimNormTbl.Date,SimNormTbl.SpotPrice,'Color',[0.9,0.9,0.9]);
hFDT = plot(SimNormTbl.Date,exp(SimNormTbl.DeterministicTrend),'y','LineWidth', 2);
hFQ = plot(SimNormTbl.Date,SimNormTbl.MCQuantiles,'Color',[0.7 0.7 0.7]);
hFM = plot(SimNormTbl.Date,SimNormTbl.MCMean,'r--');
xlabel('Date')
ylabel('Spot price ($)')
title('Electricity Spot Prices - Normal Residuals')
h = [hHSP hHDT hFSP(1) hFDT hFQ(1) hFM];
legend(h,{'Historical spot price','Historical trend','Forecasted paths',...
    'Forecasted trend','Monte Carlo quantiles','Monte Carlo mean'})
grid on

График показывает, что гораздо меньше больших скачков в предсказанных путях, выведенных от нормально распределенных невязок, чем от эпсилона, скашивает нормальные невязки.

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

simMax = max(SimTbl.SpotPrice,[],2);
simMaxNorm = max(SimNormTbl.SpotPrice,[],2);

Вычислите 6-месячные средние значения прокрутки обеих серий максимальных значений.

maxMean = movmean(simMax,183);
maxMeanNormal = movmean(simMaxNorm,183);

Постройте серию максимальных значений и соответствующих скользящих средних значений по горизонту прогноза.

figure
plot(SimTbl.Date,[simMax,simMaxNorm])
hold on
plot(SimTbl.Date,maxMean,'g','LineWidth',2)
plot(SimTbl.Date,maxMeanNormal,'m','LineWidth',2)
xlabel('Date')
ylabel('Spot price ($)')
title('Forecasted Maximum Values')
legend({'Maxima (Epsilon-Skew-Normal)','Maxima (Normal)',...
    '6-Month Moving Mean (Epsilon-Skew-Normal)',...
    '6-Month Moving Mean (Normal)'},'Location','southoutside')
grid on
hold off

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