一、奇异谱分析(Singular Spectrum Analysis, SSA) 简介
奇异谱分析(Singular Spectrum Analysis, SSA)是一种处理非线性时间序列数据的方法，通过对所要研究的时间序列的轨迹矩阵进行分解、重构等操作，提取出时间序列中的不同成分序列(长期趋势，季节趋势，噪声等)，从而进行对时间序列进行分析或去噪并用于其他一些任务。
奇异谱分析主要包括四个步骤：嵌入——分解——分组——重构。
1、嵌入，将有限长度的一维时间序列用合适的时间窗口长度L进行滞后排列得到轨迹轨迹矩阵。
2、分解，对轨迹矩阵进行奇异值分解，即分解成左矩阵、奇异值、右矩阵相乘。
3、分组，简单来说将所有的 L 个成分分组为 C 个不相交的组，代表着不同的趋势成分。
4、重构，用迟滞序列通过时间经验正交函数和时间主成分来进行重建，获得更平滑和趋势明显的序列。
详细参考：
https://blog.csdn.net/Lucky_Go/article/details/103109045
https://ssa.cf.ac.uk/zhigljavsky/pdfs/SSA/SSA_encyclopedia.pdf
二、数据准备
假设预测物品为某一股票，则在同花顺软件中下载其近二年交易数据，以收盘价作为当日价格。
同花顺数据文件格式如下：
--===============================================================
-- 文件头16个字节剖析(日线)
-- 0x6864312E3000 6 固定
-- 0x???????? 4 记录数
-- 0x4800 2 记录开始位置: 64是错的, 文件头 + 列定义 = 72 0xC000 = 192
-- 0x3800 2 每记录的长度: 56 0xB000 = 176
-- 0x0E00 2 每记录的列数: 14 0x2C00 = 44
-----------------------------------------------------------------
-- 列定义: 04表示列长度
-- 0x01300004 4 日期
-- 0x07700004 4 开盘价
-- 0x08700004 4 最高价
-- 0x09700004 4 最低价
-- 0x0B700004 4 收盘价
-- 0x13700004 4 成交金额(元)
-- 0x0D700004 4 成交量(股)
-- 0x0E700004 4 FFFFFFFF
-- 0x0F700004 4 FFFFFFFF
-- 0x11700004 4 FFFFFFFF
-- 0x12700004 4 FFFFFFFF
-- 0x50700004 4 FFFFFFFF
-- 0xE7700004 4 FFFFFFFF
-- 0xE8700004 4 FFFFFFFF
--===============================================================
其交易数据二进制文件头部数据如下：
由于其原始数据为二进制文件，需将期转换成数据表格以方便算法使用。用SQL语句将数据块直接读入数据表中。部分代码如下：
if object_id('THS.dbo.T600926') is not null drop table THS.dbo.T600926
go
declare @ varbinary(max), @max int, @e float
select @ = BulkColumn, @e = 10 from OPENROWSET(BULK N'/opt/ths/600926.day', SINGLE_BLOB) as bin
--文件记录数
select @max = substring(@,10,1)+substring(@,9,1)+substring(@,8,1)+substring(@,7,1)
--记录开始位置
select top (@max) n = identity(int,192,176) into # from syscolumns a, syscolumns b
;with zhoujy as
(
select
-- SQL没有提供按字节reverse(binary)的函数或方法，只能substring每个字节倒过来合成：
d = convert(int,substring(@, 4+n,1)+substring(@, 3+n,1)+substring(@, 2+n,1)+substring(@, 1+n,1)),
o = convert(int,substring(@, 8+n,1)+substring(@, 7+n,1)+substring(@, 6+n,1)+substring(@, 5+n,1))&0x0FFFFFFF,
p = convert(int,substring(@,12+n,1)+substring(@,11+n,1)+substring(@,10+n,1)+substring(@, 9+n,1))&0x0FFFFFFF,
q = convert(int,substring(@,16+n,1)+substring(@,15+n,1)+substring(@,14+n,1)+substring(@,13+n,1))&0x0FFFFFFF,
r = convert(int,substring(@,20+n,1)+substring(@,19+n,1)+substring(@,18+n,1)+substring(@,17+n,1))&0x0FFFFFFF,
s = convert(int,substring(@,24+n,1)+substring(@,23+n,1)+substring(@,22+n,1)+substring(@,21+n,1))&0x0FFFFFFF,
t = convert(int,substring(@,28+n,1)+substring(@,27+n,1)+substring(@,26+n,1)+substring(@,25+n,1))&0x0FFFFFFF,
u = convert(int,substring(@, 8+n,1))/16,
v = convert(int,substring(@,12+n,1))/16,
w = convert(int,substring(@,16+n,1))/16,
x = convert(int,substring(@,20+n,1))/16,
y = convert(int,substring(@,24+n,1))/16,
z = convert(int,substring(@,28+n,1))/16
from #
)
select
id = row_number()over(order by d desc),
日期 = d,
开盘 = o*power(@e,(u&7)*power(-1,sign(u&8))),
最高 = p*power(@e,(v&7)*power(-1,sign(v&8))),
最低 = q*power(@e,(w&7)*power(-1,sign(w&8))),
收盘 = r*power(@e,(x&7)*power(-1,sign(x&8))),
金额 = s*power(@e,(y&7)*power(-1,sign(y&8))),
成交 = t*power(@e,(z&7)*power(-1,sign(z&8)))
into THS.dbo.T600926 from zhoujy
得到如下数据表格式，备用
三、ML.NET中使用奇异谱算法训练模型
1、抽取数据为以下时间序列格式，以供SSA算法使用
string query = "SELECT Cast(CAST([日期] as varchar(8)) as date) as RentalDate,CAST(IIF( 日期/10000 = 2022, 1, 0 ) as Real) as Year,CAST([收盘] as REAL) as TotalRentals FROM [THS].[dbo].[T002307] order by RentalDate desc";
DatabaseSource dbSource = new DatabaseSource(SqlClientFactory.Instance,
connectionString,
query);
IDataView dataView = loader.Load(dbSource);
IDataView firstYearData = mlContext.Data.FilterRowsByColumn(dataView, "Year", upperBound: 1);
IDataView secondYearData = mlContext.Data.FilterRowsByColumn(dataView, "Year", lowerBound: 1);
从数据库读取数据，2022年 firstYearData为培训数据，2021年及以前secondYearData为验证数据。
2、定义时序分析管道
var forecastingPipeline = mlContext.Forecasting.ForecastBySsa(
outputColumnName: "ForecastedRentals",
inputColumnName: "TotalRentals",
windowSize: 7,
seriesLength: 30,
trainSize: 365,
horizon: 7,
confidenceLevel: 0.95f,
confidenceLowerBoundColumn: "LowerBoundRentals",
confidenceUpperBoundColumn: "UpperBoundRentals"); forecastingPipeline 在第一年数据中获取 365 个数据点，并按 seriesLength 参数指定的间隔从时序数据集采样或将其分为 30 天（每月）的间隔。 以一周或 7 天为一个时段分析各个样本。 确定下一个时段的预测值时，使用前面 7 天的值进行预测。 根据 horizon 参数的定义，该模型设置为预测将来的 7 个时段。 由于预测属于合理猜测，它不总是完全准确。 因此，最好了解上限和下限定义的最佳和最坏情况下的范围值。 在本案例中，设置的上下限可信度为 95%。 可信度可以相应地提高或降低。 值越高，上限和下限之间的范围越大，以便达到所需的可信度。
3、使用 Fit 方法培训模型，使数据适用于前面定义的 forecastingPipeline
SsaForecastingTransformer forecaster = forecastingPipeline.Fit(firstYearData);四、评估模型
//评做模型
Evaluate(secondYearData, forecaster, mlContext);
void Evaluate(IDataView testData, ITransformer model, MLContext mlContext)
{
IDataView predictions = model.Transform(testData);
IEnumerable<float> actual =
mlContext.Data.CreateEnumerable<ModelInput>(testData, true)
.Select(observed => observed.TotalRentals);
//获取预测值
IEnumerable<float> forecast =
mlContext.Data.CreateEnumerable<ModelOutput>(predictions, true)
.Select(prediction => prediction.ForecastedRentals[0]);
//计算实际值和预测值之间的差异，通常称为误差。
var metrics = actual.Zip(forecast, (actualValue, forecastValue) => actualValue - forecastValue);
//通过计算平均绝对误差和均方根误差值来衡量性能。
var MAE = metrics.Average(error => Math.Abs(error)); // Mean Absolute Error
var RMSE = Math.Sqrt(metrics.Average(error => Math.Pow(error, 2))); // Root Mean Squared Error
//将指标输出到控制台。
Console.WriteLine("Evaluation Metrics");
Console.WriteLine("---------------------");
Console.WriteLine($"Mean Absolute Error: {MAE:F3}");
Console.WriteLine($"Root Mean Squared Error: {RMSE:F3}\n");
}
保存模型供项目使用
//保存模型
var forecastEngine = forecaster.CreateTimeSeriesEngine<ModelInput, ModelOutput>(mlContext);
forecastEngine.CheckPoint(mlContext, modelPath);
五、使用模型预测结果
//使用模型
Forecast(secondYearData, 7, forecastEngine, mlContext);
void Forecast(IDataView testData, int horizon, TimeSeriesPredictionEngine<ModelInput, ModelOutput> forecaster, MLContext mlContext)
{
ModelOutput forecast = forecaster.Predict();
IEnumerable<string> forecastOutput =
mlContext.Data.CreateEnumerable<ModelInput>(testData, reuseRowObject: false)
.Take(horizon)
.Select((ModelInput rental, int index) =>
{
string rentalDate = rental.RentalDate.ToShortDateString();
float actualRentals = rental.TotalRentals;
float lowerEstimate = Math.Max(0, forecast.LowerBoundRentals[index]);
float estimate = forecast.ForecastedRentals[index];
float upperEstimate = forecast.UpperBoundRentals[index];
return $"Date: {rentalDate}\n" +
$"Actual Rentals: {actualRentals}\n" +
$"Lower Estimate: {lowerEstimate}\n" +
$"Forecast: {estimate}\n" +
$"Upper Estimate: {upperEstimate}\n";
});
// Output predictions
Console.WriteLine("Rental Forecast");
Console.WriteLine("---------------------");
foreach (var prediction in forecastOutput)
{
Console.WriteLine(prediction);
}
}
运行结果
六、参考文档
1、Tutorial: Forecast bike rental demand - time series - ML.NET | Microsoft Learn
2、https://blog.csdn.net/Lucky_Go/article/details/103109045
3、https://ssa.cf.ac.uk/zhigljavsky/pdfs/SSA/SSA_encyclopedia.pdf
