对于命令+插值：需要一些技巧

Question

我有一个矩阵A，其中包含三列：每日日期、价格和小时数 - 都是相同大小的向量 - 一天中有多个价格与小时数相关联。

下面的示例数据：

A_dates =          A_hours=        A_prices=
[20080902         [9.698           [24.09
20080902          9.891             24.59
200080902         10.251            24.60 
20080903          9.584             25.63
200080903         10.45             24.96
200080903         12.12             24.78
200080904          12.95            26.98 
20080904           13.569           26.78
20080904]          14.589]          25.41]

• 请记住，我有大约两年的每日数据，每天大约有 10000 个价格，几乎涵盖了每天上午 9:30 到下午 16:00 的每一分钟。实际上我最初的数据集时间是以毫秒为单位的。然后我以小时为单位转换我的毫秒数。我有几个小时，例如 14.589，以 3 种不同的价格重复了 3 次。因此我做了以下事情：

  时间=[A_dates,A_hours,A_prices]; [timeinhr,price]=consolidator(time,A_prices,'mean');其中 timeinhr 是向量 A_dates 和 A_hours

取平均价格，例如 14.589 小时。 然后对于任何缺少 .25 .50 .75 和整数小时的小时数 - 我希望插值。

对于每个日期，小时重复，我需要线性插入一些“想要”小时没有的价格。但是，如果我的时间在我的专栏中重复，我当然不能使用命令 interp1，因为我有好几天。所以说：

%# here I want hours in 0.25unit increments (like 9.5hrs)
new_timeinhr = 0:0.25:max(A_hours));

day_hour = rem(new_timeinhour, 24);

%# Here I want only prices between 9.5hours and 16hours
new_timeinhr( day_hour <= 9.2 | day_hour >= 16.1 ) = [];  

然后我创建了一个唯一的一天向量，并希望使用 for 和 if 命令插入 daily，然后将我的新价格一个接一个地堆叠在一个向量中：

days = unique(A_dates);
for j = 1:length(days);
    if A_dates == days(j)
       int_prices(j) = interp1(A_hours, A_prices, new_timeinhr);
    end;
end;

我的错误是：

In an assignment A(I) = B, the number of elements in B and I must be the same.

如何将int_prices(j) 写入堆栈？

Answer 1

我建议将您的输入转换为单个单调时间值。使用 MATLAB datenum 格式，将一天表示为 1。这样做有很多优点：您可以获得内置的 MATLAB 时间/日期函数，您可以通过 datetick 获得格式化为日期/时间的绘图标签，以及插值只是工作。没有测试数据，我无法测试这段代码，但这是大致的思路。

根据您将日期存储为 20080902（我假设为 yyyymmdd）的新信息，我更新了初始转换代码。另外，由于 A 的布局会引起混乱，我将把 A 的列称为向量 A_prices、A_hours 和 A_dates。

% This datenum vector matches A.  I'm assuming they're already sorted by date and time
At = datenum(num2str(A_dates), 'yyyymmdd') + datenum(0, 0, 0, A_hours, 0, 0);
incr = datenum(0, 0, 0, 0.25, 0, 0);  % 0.25 hour
t = (At(1):incr:At(end)).';       % Full timespan of dataset, in 0.25 hour increments

frac_hours = 24*(t - floor(t));        % Fractional hours into the day
t_business_day = t((frac_hours > 9.4) & (frac_hours < 16.1));  % Time vector only where you want it

P = interp1(At, A_prices, t_business_day);

我再说一遍，因为没有测试数据，我无法测试代码。我强烈建议使用datestr 将日期转换代码back 从 datenum 转换为可读日期来测试日期转换代码。

Answer 2

按照@Peter 的建议，将天数/小时数转换为序列日期数字绝对是可行的方法。根据他的代码（我已经赞成），我在下面给出一个简单的例子。

首先，我创建一些类似于您所描述的假数据（也有一些缺失的部分）：

%# three days in increments of 1 hour
dt = datenum(num2str((0:23)','2012-06-01 %02d:00'), 'yyyy-mm-dd HH:MM');   %#'
dt = [dt; dt+1; dt+2];

%# price data corresponding to each hour
p = cumsum(rand(size(dt))-0.5);

%# show plot
plot(dt, p, '.-'), datetick('x')
grid on, xlabel('Date/Time'), ylabel('Prices')

%# lets remove some rows as missing
idx = ( rand(size(dt)) < 0.1 );
hold on, plot(dt(idx), p(idx), 'ro'), hold off
legend({'prices','missing'})
dt(idx) = [];
p(idx) = [];

%# matrix same as yours: days,prices,hours
ymd = str2double( cellstr(datestr(dt,'yyyymmdd')) );
hr = str2double( cellstr(datestr(dt,'HH')) );
A = [ymd p hr];

%# let clear all variables except the data matrix A
clearvars -except A

接下来我们以 15 分钟为增量在整个范围内插入价格数据：

%# convert days/hours to serial date number
dt = datenum(num2str(A(:,[1 3]),'%d %d'), 'yyyymmdd HH');

%# create a vector of 15 min increments
t_15min = (0:0.25:(24-0.25))';                  %#'
tt = datenum(0,0,0, t_15min,0,0);

%# offset serial date across all days
ymd = datenum(num2str(unique(A(:,1))), 'yyyymmdd');
tt = bsxfun(@plus, ymd', tt);                   %#'
tt = tt(:);

%# interpolate data at new datetimes
pp = interp1(dt, A(:,2), tt);

%# extract desired period of time from each day
idx = (9.5 <= t_15min & t_15min <= 16);
idx2 = bsxfun(@plus, find(idx), (0:numel(ymd)-1)*numel(t_15min));
P = pp(idx2(:));

%# plot interpolated data, and show extracted periods
figure, plot(tt, pp, '.-'), datetick('x'), hold on
plot([tt(idx2);nan(1,numel(ymd))], [pp(idx2);nan(1,numel(ymd))], 'r.-')
hold off, grid on, xlabel('Date/Time'), ylabel('Prices')
legend({'interpolated prices','period of 9:30 - 16:00'})

以下是显示原始数据和插值数据的两个图：

Answer 3

我想我可以这样解决它：

new_timeinhr = 0:0.25:max(A(:,2));
day_hour = rem(new_timeinhr, 24);
new_timeinhr( day_hour <= 9.4 | day_hour >= 16.1 ) = [];

days=unique(data(:,1));
P=[];
for j=1:length(days);
    condition=A(:,1)==days(j);
    intprices = interp1(A(condition,2), A(condition,3), new_timeinhr);
    P=vertcat(P,intprices');
end;