Обзор модели и анализа

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

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

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

Загрузка и визуализация спотовых цен на электроэнергию

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

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

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

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

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-месячному циклам.

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

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

где:

Синусоидальные и косинусные члены в модели учитывают возможное фазовое смещение. Это означает, что для фазового смещения$$$$

$$\mathrm{Asin}(ft+start)\mathrm{=}(\mathrm{}\mathrm{Acosü})\mathrm{sinft}+\mathrm{(}\mathrm{Asinstart}\mathrm{)}\mathrm{cosft}\mathrm{=}$$Bsinft + Ccosft,

где $$B=\mathrm{}$$Acos $$\mathrm{}$$Поэтому включение синусоидальных и косинусных членов с одной и той же частотой является причиной фазового сдвига.

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

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 сбрасывает $$\mathrm{член\; sin}\mathrm{(}2xeont$$) с линейной модели, указывая, что годовой цикл начинается с волнового пика.

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

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}\mathrm{=}$0. Эти особенности предполагают, что авторегрессионная модель подходит для отклоненных данных.

Чтобы определить количество лагов для включения в авторегрессионную модель, примените методологию Бокса-Дженкинса. Постройте график автокорреляционной функции (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) с формой

$${}_{\mathrm{startt}}={}_{\mathrm{}}{\mathrm{startt-1}}_{}$$+ ut,

где $$\mathrm{start}$$- авторегрессионный коэффициент запаздывания 1, а $${}_{}\mathrm{ut\sim N}\mathrm{(}0{,}^{\mathrm{}}$$start2) - 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-на-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

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

Серия отклоненных остатков модели

Распределение эпсилон-перекос-нормаль представляет собой семейство распределения, близкое к нормальному, с расположением, $$шкалой, $и$дополнительным параметром перекосов $$Параметр skewness моделирует любой ненулевой перекос в данных [2]. Если λ $$=0$$, распределение эпсилон-перекос-нормаль уменьшается до нормального распределения.

Эконометрика 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')])

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

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

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

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

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

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

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

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

Смоделировать 1000 путей сниженных цен журнала путем фильтрации смоделированных стандартизированных остатков через оценочную модель AR (1)EstDTMdl. Определение расчетных условных отклонений condVar в качестве предварительного примера условных отклонений наблюдаемые сниженные цены журнала в DataTable в качестве предварительных ответов и оценочных стандартизированных остатков в DataTable в качестве предварительных нарушений. В этом контексте «presample» - это период до горизонта прогноза.

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 - долгосрочный детерминированный трендовый компонент спотовых цен журнала в горизонте прогноза.

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

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,5-го, 5-го, 25-го, 75-го, 95-го и 97,5-го процентилей.

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

Кроме того, при наличии лицензии на финансовый Toolbox™ можно использовать fanplot для построения графика квантилей Монте-Карло.

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

Рассчитайте ежемесячные значения Monte Carlo и стандартные отклонения моделируемых путей спотовых цен.

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

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

Сравнение результатов Epsilon-Skew-Normal с предположением о нормальном распределении

В предыдущем анализе стандартизированные остатки получаются из установленного распределения эпсилон-перекос-нормаль. Получение смоделированных спотовых цен путем обратного преобразования нормально распределенных стандартизированных остатков. Стандартизированные остатки можно моделировать непосредственно с помощью 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,5-го, 5-го, 25-го, 75-го, 95-го и 97,5-го процентилей.

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

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