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

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

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

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

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

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

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

Рабочая область содержит расписание MATLAB ® 1411 на 1 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\cos \theta )\sin ft+(A\sin \theta )\cos ft=B\sin ft+C\cos ft$$,

где $$B=A\cos \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$. Эти функции предполагают, что авторегрессионная модель подходит для детрендированных данных.

Чтобы определить количество лагов, которые будут включены в авторегрессивную модель, примените методологию Box-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, an 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$. Параметр skewness моделирует любой ненулевой перекос в данных [2]. Если $$\epsilon =0$$, распределение эпсилон-наклон-нормаль уменьшается до нормального распределения.

Econometrics Toolbox™ включает пользовательский объект распределения, представляющий распределение epsilon-skew-normal в 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 установленного распределения epsilon-skew-normal. Верните точки, в которых программное обеспечение оценивает эмпирический 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')])

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

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

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

Сохраните данные моделирования в расписании с именем 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-by- numPaths матрица моделируемых стандартизированных остаточных путей. Строки соответствуют периодам в прогнозном горизонте, а столбцы соответствуют отдельным путям симуляции.

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

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

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

Кроме того, если у вас есть лицензия Financial 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

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

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

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

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