Модель и симулирует спотовые цены электричества Используя Скошенное Нормальное распределение

В этом примере показано, как симулировать будущее поведение спотовых цен электричества из модели временных рядов, подбиравшей к историческим данным. Поскольку спотовые цены электричества могут показать большие отклонения, пример моделирует инновации с помощью скошенного нормального распределения. Набор данных в этом примере содержит симулированные, ежедневные спотовые цены электричества с 1 января 2010 до 11 ноября 2013.

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

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

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

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

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

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

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

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

  5. Задайте и оцените соответствующую авторегрессивную модель (ARIMA) скользящего среднего значения для детрендированных данных путем применения методологии Поля-Jenkins.

  6. Соответствуйте скошенному нормальному распределению вероятностей к стандартизированным остаточным значениям подбиравшей модели ARIMA. Этот шаг требует, чтобы пользовательский объект вероятностного распределения создал использование среды, доступной в Statistics and Machine Learning Toolbox™.

  7. Симулируйте спотовые цены. Во-первых, чертите iid случайный стандартизированный остаточный ряд от строившего распределения вероятности от шага 6. Затем backtransform симулированные остаточные значения путем инвертирования шагов 1-5.

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

Загрузите 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

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

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

  1. Регресс логарифмическая спотовая цена на номер лет протек начиная с запуска данных при помощи polyfit.

  2. Оцените среднюю логарифмическую спотовую цену для каждой временной стоимости при помощи polyval.

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'})

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

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

  2. Применяйте преобразование Фурье к результату.

  3. Вычислите спектр мощности преобразованных данных и соответствующего вектора частоты.

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

Подбирайте долгосрочную детерминированную модель тренда

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

LogPrice=β0+β1t+β2sin(2πt)+β3cos(2πt)+β4sin(4πt)+β5cos(4πt)+ξt,

где:

  • t годы, истекшие от запуска выборки.

  • ξtN(0,τ2) iid последовательность случайных переменных.

  • β2 и β3 коэффициенты ежегодных компонентов.

  • β4 и β5 коэффициенты полугодовых компонентов.

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

Asin(ft+θ)=(Acosθ)sinft+(Asinθ)cosft=Bsinft+Ccosft,

где B=Acosθ и C=Asinθ константы. Поэтому включение синуса и условий косинуса с той же частотой составляет смещение фазы.

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

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π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

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

ξt=ϕξt-1+ut,

где ϕ задержка 1 авторегрессивный коэффициент и utN(0,σ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

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

Скошенный остаточный ряд модели

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

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

  1. Откройте приложение Distribution Fitter (Statistics and Machine Learning Toolbox™) путем ввода distributionFitter в командной строке.

  2. В приложении выберите File> Define Custom Distributions. Редактор отображает файл определения класса, содержащий выборку класса пользовательского дистрибутива (распределение Лапласа). Закройте окно Define Custom Distributions и окно Distribution Fitter.

  3. В демонстрационном файле определения класса измените имя класса от LaplaceDistribution к EpsilonSkewNormal.

  4. В вашей текущей папке создайте директорию под названием +prob, затем сохраните новый файл определения класса в той папке.

  5. В файле определения класса введите параметры скошенного нормального распределения эпсилона как свойства и настройте методы так, чтобы класс представлял скошенное нормальное распределение эпсилона. Для получения дополнительной информации на распределении, см. [2].

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

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]

Предполагаемый параметр скошенности εˆ приблизительно-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   

εˆ имеет небольшую стандартную погрешность (~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

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

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

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

  1. Сравните эмпирическое, и соответствовал (ссылочным) кумулятивным функциям распределения (cdfs).

  2. Проведите тест Кольмогорова-Смирнова для качества подгонки.

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

Вычислите эмпирический 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 кажутся подобными.

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

  1. Если данные, которые включают эмпирический cdf, являются тем же набором как данные, которые включают ссылку cdf, то тест Кольмогорова-Смирнова недопустим. Поэтому разделите набор данных в наборы обучающих данных и наборы тестов при помощи инструментов раздела перекрестной проверки Statistics and Machine Learning Toolbox™.

  2. Соответствуйте распределению к набору обучающих данных.

  3. Подтвердите подгонку на наборе тестов.

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

  1. Задайте даты горизонта прогноза.

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

  3. Получите симулированные, детрендированные цены журнала путем пропущения остаточных значений через подбиравшую модель AR (1).

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

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

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

Сохраните данные моделирования в расписании под названием 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 построить квантили Монте-Карло.

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

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

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

  3. Постройте график лучшей подгонки к объединенному историческому и моменты Монте-Карло.

  4. Визуально сравните линии лучшей подгонки.

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

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

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

Ссылки

[1] Люсия, J.J., и Э. Шварц. "Цены на электроэнергию и производные степени: Доказательство от скандинавского Exchange Степени". Анализ Исследования Производных. Издание 5, № 1, 2002, стр 5–50.

[2] Mudholkara, G.S., и нашей эры Хатсон. "Скошенное Нормальное распределение Эпсилона для Анализа Почти нормальных Данных". Журнал Статистического Планирования и Вывода. Издание 83, № 2, 2000, стр 291–309.

Смотрите также

Объекты

Функции

Похожие темы

Для просмотра документации необходимо авторизоваться на сайте