В этом примере показано, как симулировать будущее поведение спотовых цен электричества из модели временных рядов, подбиравшей к историческим данным. Поскольку спотовые цены электричества могут показать большие отклонения, пример демонстрирует инновации с помощью скошенного нормального распределения. Набор данных в этом примере содержит симулированные, ежедневные спотовые цены электричества с 1 января 2010 до 11 ноября 2013.
Цена электричества сопоставлена с соответствующим спросом [1]. Изменения в численности населения и техническом прогрессе предлагают, чтобы спотовые цены электричества имели долгосрочный тренд. Кроме того, учитывая климат области, спрос на электричество колеблется согласно сезону. Следовательно, спотовые цены электричества колеблются, точно так же предлагая включение сезонного компонента в модели.
В периоды высокого спроса спотовые цены могут показать скачки, когда технический персонал добавляет текущее предоставление с энергией, произведенной от менее эффективных методов. Эти периоды высокого спроса предполагают, что инновационное распределение спотовых цен электричества правильно скошенный, а не симметричный.
Характеристики спотовых цен электричества влияют на шаги для создания модели:
Определите, имеет ли ряд спотовой цены экспоненциальные, долгосрочные детерминированные, и сезонные тренды путем осмотра его графика временных рядов и выполнения других статистических тестов.
Если ряд спотовой цены показывает экспоненциальный тренд, удалите его путем применяния логарифмического преобразования к данным.
Оцените, имеют ли данные долгосрочный детерминированный тренд. Идентифицируйте линейную линейную регрессию использования компонента. Идентифицируйте доминирующие сезонные компоненты путем применения спектрального анализа к данным.
Выполните пошаговую линейную регрессию, чтобы оценить полный долгосрочный детерминированный тренд. Получившийся тренд включает линейные и сезонные компоненты с частотами, определенными спектральным анализом. Детрендируйте ряд путем удаления предполагаемого долгосрочного детерминированного тренда из данных.
Задайте и оцените соответствующую авторегрессивную модель (ARIMA) скользящего среднего значения для детрендированных данных путем применения методологии Поля-Jenkins.
Соответствуйте скошенному нормальному распределению вероятностей к стандартизированным остаточным значениям подбиравшей модели ARIMA. Этот шаг требует, чтобы пользовательский объект вероятностного распределения создал использование среды, доступной в Statistics and Machine Learning Toolbox™.
Симулируйте спотовые цены. Во-первых, чертите 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
Поскольку график спотовых цен клонится вниз и имеет сезонный шаблон, долгосрочный тренд может содержать линейные и сезонные компоненты. Оценить линейный компонент тренда:
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-месячным циклам.
Путем объединения доминирующих сезонных компонентов, оцененных спектральным анализом, с линейным компонентом, оцененным наименьшими квадратами, линейная модель для логарифмических спотовых цен может иметь эту форму:
где:
годы, истекшие от запуска выборки.
iid последовательность случайных переменных.
и коэффициенты ежегодных компонентов.
и коэффициенты полугодовых компонентов.
Синус и члены в модели косинуса составляют возможное смещение фазы. Таким образом, поскольку фаза возмещена :
,
где и константы. Поэтому включение синуса и условий косинуса с той же частотой составляет смещение фазы.
Сформируйте матрицу проекта. Поскольку программное обеспечение включает постоянный член в модели по умолчанию, не включайте столбец из единиц в матрице проекта.
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
понижается назовите из линейной модели, указав, что ежегодный цикл начинается на пике волны.
Постройте подобранную модель.
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
Детрендированные данные кажутся в центре в нуле, и ряд показывает последовательную автокорреляцию, потому что несколько кластеров последовательных остаточных значений происходят выше и ниже . Эти функции предполагают, что авторегрессивная модель подходит для детрендированных данных.
Чтобы определить количество задержек, чтобы включать в авторегрессивную модель, примените методологию Поля-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) процесс с формой
,
где задержка 1 авторегрессивный коэффициент и 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]. Если , скошенное нормальное распределение эпсилона уменьшает до нормального распределения.
Econometrics Toolbox™ включает объект пользовательского дистрибутива представление скошенного нормального распределения эпсилона в mlr/examples/econ/helper/+prob/EpsilonSkewNormalDistribution.m
, где mlr
значение matlabroot
. Чтобы создать ваш собственный объект распределения, выполните эти шаги:
Откройте приложение Distribution Fitter (Statistics and Machine Learning Toolbox™) путем ввода distributionFitter
в командной строке.
В приложении выберите File> Define Custom Distributions. Редактор отображает файл определения класса, содержащий выборку класса пользовательского дистрибутива (распределение Лапласа). Закройте окно Define Custom Distributions и окно Distribution Fitter.
В демонстрационном файле определения класса измените имя класса от LaplaceDistribution
к EpsilonSkewNormal
.
В вашей текущей папке создайте директорию под названием +prob
, затем сохраните новый файл определения класса в той папке.
В файле определения класса введите параметры скошенного нормального распределения эпсилона как свойства и настройте методы так, чтобы класс представлял скошенное нормальное распределение эпсилона. Для получения дополнительной информации на распределении, см. [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
Подходящее распределение, кажется, хорошая модель для остаточных значений.
Оцените качество подгонки скошенного нормального распределения эпсилона путем выполнения этой процедуры:
Сравните эмпирическое, и соответствовал (ссылочным) кумулятивным функциям распределения (cdfs).
Проведите тест Кольмогорова-Смирнова для качества подгонки.
Постройте максимум 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 кажутся подобными.
Проведите тест Кольмогорова-Смирнова для качества подгонки путем выполнения этой процедуры:
Если данные, которые включают эмпирический cdf, являются тем же набором как данные, которые включают ссылку cdf, то тест Кольмогорова-Смирнова недопустим. Поэтому разделите набор данных в наборы обучающих данных и наборы тестов при помощи инструментов раздела перекрестной проверки Statistics and Machine Learning Toolbox™.
Соответствуйте распределению к набору обучающих данных.
Подтвердите подгонку на наборе тестов.
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 года после конца исторических данных. Создайте модель, из которой можно симулировать, состоявший из предполагаемых компонентов временных рядов, путем выполнения этой процедуры:
Задайте даты горизонта прогноза.
Получите симулированные остаточные значения путем симуляции стандартизированных остаточных значений подходящего скошенного нормального распределения эпсилона, затем масштабирования результата предполагаемыми мгновенными стандартными отклонениями.
Получите симулированные, детрендированные цены журнала путем пропущения остаточных значений через подбиравшую модель AR (1).
Предскажите значения детерминированного тренда с помощью подобранной модели.
Получите симулированные цены журнала путем объединения симулированных, детрендированных цен журнала и предсказанных, детерминированных значений тренда.
Получите симулированные спотовые цены возведением в степень симулированные логарифмические спотовые цены.
Сохраните данные моделирования в расписании под названием 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
Несмотря на то, что серии симулированных максимумов коррелируются, график показывает, что эпсилон скашивается, нормальные максимумы больше, чем нормальные максимумы.
[1] Люсия, J.J., и Э. Шварц. "Цены на электроэнергию и производные степени: Доказательство от скандинавского Exchange Степени". Анализ Исследования Производных. Издание 5, № 1, 2002, стр 5–50.
[2] Mudholkara, G.S., и нашей эры Хатсон. "Скошенное Нормальное распределение Эпсилона для Анализа Почти нормальных Данных". Журнал Статистического Планирования и Вывода. Издание 83, № 2, 2000, стр 291–309.