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

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

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

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

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

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

Загрузите Data_ElectricityPrices MAT-файл включен с 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- 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- 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

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