精华内容
下载资源
问答
  • 目录1 数据相关2 时间序列中的模型(Patterns)3 如何分解时间序列中的各个成分4 平稳与不平稳时间序列4.1 这些数据有什么明显的特点?4.2 为什么要在预测前把序列变成平稳的?4.3 如何对平稳性进行测试4.4 白噪音和...


    1 数据相关


    2 时间序列中的模型(Patterns)

    任何时间序列都能够被分解成下面几个部分:
    BaseLevel+Trend+Seasonality+Error Base Level + Trend + Seasonality + Error
    Additive time series:
    Value=BaseLevel+Trend+Seasonality+Error Value = Base Level + Trend + Seasonality + Error
    Multiplicative Time Series:
    Value=BaseLevelTrendSeasonalityError Value = Base Level * Trend * Seasonality * Error


    3 如何分解时间序列中的各个成分

    我们可以将一个时间序列看成是由base level,trend,seasonal index,the residual 这些成分组成的Additive time seriesMultiplicative Time Series,从而使用传统的分解方法。

    我们可以使用statsmodels 库里的seasonal_decompose() 函数来进行分解。

    from statsmodels.tsa.seasonal import seasonal_decompose
    from dateutil.parser import parse
    
    # Import Data
    df = pd.read_csv('https://raw.githubusercontent.com/selva86/datasets/master/a10.csv', parse_dates=['date'], index_col='date')
    
    # Multiplicative Decomposition 
    result_mul = seasonal_decompose(df['value'], model='multiplicative', extrapolate_trend='freq')
    
    # Additive Decomposition
    result_add = seasonal_decompose(df['value'], model='additive', extrapolate_trend='freq')
    
    # Plot
    plt.rcParams.update({'figure.figsize': (10,10)})
    result_mul.plot().suptitle('Multiplicative Decompose', fontsize=22)
    result_add.plot().suptitle('Additive Decompose', fontsize=22)
    plt.show()
    

    分解因素

    # Extract the Components ----
    # Actual Values = Product of (Seasonal * Trend * Resid)
    df_reconstructed = pd.concat([result_mul.seasonal, result_mul.trend, result_mul.resid, result_mul.observed], axis=1)
    df_reconstructed.columns = ['seas', 'trend', 'resid', 'actual_values']
    df_reconstructed.head()
    
    [output]:
    ADF Statistic: 3.1451856893067434
    p-value: 1.0
    Critial Values:
       1%, -3.465620397124192
    Critial Values:
       5%, -2.8770397560752436
    Critial Values:
       10%, -2.5750324547306476
    
    KPSS Statistic: 1.313675
    p-value: 0.010000
    Critial Values:
       10%, 0.347
    Critial Values:
       5%, 0.463
    Critial Values:
       2.5%, 0.574
    Critial Values:
       1%, 0.739
    

    4 平稳与不平稳时间序列

    4.1 这些数据有什么明显的特点?

    在这里插入图片描述

    4.2 为什么要在预测前把序列变成平稳的?

    最简单的原因是:平稳序列的预测更简单,而且预测结果更可靠。

    一个重要原因是:自回归预测模型是一种线性模型,它会用到序列自身作为预测因子(Predictors),也就是自变量。我们知道,线性回归当自变量彼此之间没有相关关系的时候,表现最好。因此,我们需要让序列中的每个因素相互独立。

    4.3 如何对平稳性进行测试

    第一种方法:可视化。平稳序列就像之前的图片中展示的那样;非平稳序列同样。

    第二种方法:将序列且分成2 或多个相邻的序列,然后计算每个序列的统计量,比如均值、方差,自相关系数。如果每段的这些统计量差别很大,那么这个序列很可能是非平稳的。

    第三种方法:ADH Test

    4.4 白噪音和平稳序列的区别

    • 相同:关于时间的函数,均值和方差都不会随时间变化。
    • 不同:关于0 均值的完全随机的波动。

    5 如何去掉时间序列中的趋势

    1. 用观测数据减去拟合的趋势线。拟合的趋势线可能包含一个线性回归的模型,该模型的一个自变量是时间戳数据。更复杂的趋势,我们可以考虑加入二次项。
    2. 第2节中,用观测数据减去分解的时间序列中的趋势部分。
    3. 减去均值。
    4. 应用滤波器。比如Baxter-King filter(statsmodels.tsa.filters.bkfilter)

    前两种方法如下所示。

    # Using scipy: Subtract the line of best fit
    from scipy import signal
    df = pd.read_csv('https://raw.githubusercontent.com/selva86/datasets/master/a10.csv', parse_dates=['date'])
    detrended = signal.detrend(df.value.values)
    plt.plot(detrended)
    plt.title('Drug Sales detrended by subtracting the least squares fit', fontsize=16)
    

    detrend

    # Using statmodels: Subtracting the Trend Component.
    from statsmodels.tsa.seasonal import seasonal_decompose
    df = pd.read_csv('https://raw.githubusercontent.com/selva86/datasets/master/a10.csv', parse_dates=['date'], index_col='date')
    result_mul = seasonal_decompose(df['value'], model='multiplicative', extrapolate_trend='freq')
    detrended = df.value.values - result_mul.trend
    plt.plot(detrended)
    plt.title('Drug Sales detrended by subtracting the trend component', fontsize=16)
    

    detrend


    6 如何去掉时间序列中的季节

    1. 用和季节长度一致的窗口进行滑动平均。
    2. 用观测数据减去前一季节数据。
    3. Divide the series by the seasonal index obtained from STL decomposition
    # Subtracting the Trend Component.
    df = pd.read_csv('https://raw.githubusercontent.com/selva86/datasets/master/a10.csv', parse_dates=['date'], index_col='date')
    
    # Time Series Decomposition
    result_mul = seasonal_decompose(df['value'], model='multiplicative', extrapolate_trend='freq')
    
    # Deseasonalize
    deseasonalized = df.value.values / result_mul.seasonal
    
    # Plot
    plt.plot(deseasonalized)
    plt.title('Drug Sales Deseasonalized', fontsize=16)
    plt.plot()
    

    deseasonality述

    6.1 怎么测试序列中的季节性?

    如果我们想要对季节性更精准的检测,那么就是用Autocorrelation Function (ACF) plot。当ACF 中有很强的季节性模式,ACF plot 通常会表现出明显的按照季节规律的峰值

    比如在我们的数据(Drug Sales)中,我们能够发现每年的季节性,12th,24th,36th……


    7 缺失值的处理

    如果数据测量值本身为0,这种情况只需要填充0。

    如果是数据缺失的情况,对于平稳序列,那么一般可以填充均值;对于非平稳序列,我们不应该用均值,比较简单粗暴的办法是填充前一个观测值。

    有几种比较有效的办法:

    • 向后填充法
    • 线性插值法??
    • 二次插值法??
    • 最近邻均值法
    • 季节均值法

    为了评估缺失值的填充效果,我在时间序列中手动加入缺失值,用以上几种方法对其进行填充,然后计算填充后的序列与原序列的均方误差。

    # # Generate dataset
    from scipy.interpolate import interp1d
    from sklearn.metrics import mean_squared_error
    df_orig = pd.read_csv(r'D:\jupyter files\data_practice_python\a10.csv', parse_dates=['date'], index_col='date')
    df = pd.read_csv(r'D:\jupyter files\data_practice_python\a10_missings.csv', parse_dates=['date'], index_col='date')
    
    fig, axes = plt.subplots(7, 1, sharex=True, figsize=(10, 12))
    plt.rcParams.update({'xtick.bottom' : False})
    
    ## 1. Actual -------------------------------
    df_orig.plot(title='Actual', ax=axes[0], label='Actual', color='red', style=".-")
    df.plot(title='Actual', ax=axes[0], label='Actual', color='green', style=".-")
    axes[0].legend(["Missing Data", "Available Data"])
    
    ## 2. Forward Fill --------------------------
    df_ffill = df.ffill()
    error = np.round(mean_squared_error(df_orig['value'], df_ffill['value']), 2)
    df_ffill['value'].plot(title='Forward Fill (MSE: ' + str(error) +")", ax=axes[1], label='Forward Fill', style=".-")
    
    ## 3. Backward Fill -------------------------
    df_bfill = df.bfill()
    error = np.round(mean_squared_error(df_orig['value'], df_bfill['value']), 2)
    df_bfill['value'].plot(title="Backward Fill (MSE: " + str(error) +")", ax=axes[2], label='Back Fill', color='firebrick', style=".-")
    
    ## 4. Linear Interpolation ------------------
    df['rownum'] = np.arange(df.shape[0])
    df_nona = df.dropna(subset = ['value'])
    f = interp1d(df_nona['rownum'], df_nona['value'])
    df['linear_fill'] = f(df['rownum'])
    error = np.round(mean_squared_error(df_orig['value'], df['linear_fill']), 2)
    df['linear_fill'].plot(title="Linear Fill (MSE: " + str(error) +")", ax=axes[3], label='Cubic Fill', color='brown', style=".-")
    
    ## 5. Cubic Interpolation --------------------
    f2 = interp1d(df_nona['rownum'], df_nona['value'], kind='cubic')
    df['cubic_fill'] = f2(df['rownum'])
    error = np.round(mean_squared_error(df_orig['value'], df['cubic_fill']), 2)
    df['cubic_fill'].plot(title="Cubic Fill (MSE: " + str(error) +")", ax=axes[4], label='Cubic Fill', color='red', style=".-")
    
    # Interpolation References:
    # https://docs.scipy.org/doc/scipy/reference/tutorial/interpolate.html
    # https://docs.scipy.org/doc/scipy/reference/interpolate.html
    
    # ## 6. Mean of 'n' Nearest Past Neighbors ------
    # def knn_mean(ts, n):
    #     out = np.copy(ts)
    #     for i, val in enumerate(ts):
    #         if np.isnan(val):
    #             n_by_2 = np.ceil(n/2)
    #             lower = np.max([0, int(i-n_by_2)])
    #             upper = np.min([len(ts)+1, int(i+n_by_2)])
    #             ts_near = np.concatenate([ts[lower:i], ts[i:upper]])
    #             out[i] = np.nanmean(ts_near)
    #     return out
    
    # df['knn_mean'] = knn_mean(df.value.values, 8)
    # error = np.round(mean_squared_error(df_orig['value'], df['knn_mean']), 2)
    # df['knn_mean'].plot(title="KNN Mean (MSE: " + str(error) +")", ax=axes[5], label='KNN Mean', color='tomato', alpha=0.5, style=".-")
    
    ## 7. Seasonal Mean ----------------------------
    def seasonal_mean(ts, n, lr=0.7):
        """
        Compute the mean of corresponding seasonal periods
        ts: 1D array-like of the time series
        n: Seasonal window length of the time series
        """
        out = np.copy(ts)
        for i, val in enumerate(ts):
            if np.isnan(val):
                ts_seas = ts[i-1::-n]  # previous seasons only
                if np.isnan(np.nanmean(ts_seas)):
                    ts_seas = np.concatenate([ts[i-1::-n], ts[i::n]])  # previous and forward
                out[i] = np.nanmean(ts_seas) * lr
        return out
    
    df['seasonal_mean'] = seasonal_mean(df.value, n=12, lr=1.25)
    error = np.round(mean_squared_error(df_orig['value'], df['seasonal_mean']), 2)
    df['seasonal_mean'].plot(title="Seasonal Mean (MSE: " + str(error) +")", ax=axes[6], label='Seasonal Mean', color='blue', alpha=0.5, style=".-")
    

    8 自相关和偏自相关函数

    8.1 概念

    详见 https://blog.csdn.net/qq_41103204/article/details/105810742

    from statsmodels.tsa.stattools import acf, pacf
    from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
    
    df = pd.read_csv('https://raw.githubusercontent.com/selva86/datasets/master/a10.csv')
    
    # Calculate ACF and PACF upto 50 lags
    # acf_50 = acf(df.value, nlags=50)
    # pacf_50 = pacf(df.value, nlags=50)
    
    # Draw Plot
    fig, axes = plt.subplots(1,2,figsize=(16,3), dpi= 100)
    plot_acf(df.value.tolist(), lags=50, ax=axes[0])
    plot_pacf(df.value.tolist(), lags=50, ax=axes[1])
    

    在这里插入图片描述

    8.2 如何计算偏自相关系数?

    序列滞后kk 处的偏自相关系数是yy 的自回归方程的滞后系数。yy 的自回归方程是指 yy 以自己的滞后作为预测因子的线性回归。

    举个例子,如果yty_{t} 为当前序列,yt1y_{t-1} 即为滞后期为1 的yy 值,那么滞后期为 3 处yt3y_{t-3} 的偏自相关系数是下面方程中的 α3\alpha_{3}

    yt=α0+α1yt1+α2yt2+α3yt3 y_{t} = \alpha_{0} + \alpha_{1} y_{t-1} + \alpha_{2} y_{t-2} + \alpha_{3} y_{t-3}

    8.3 滞后图

    滞后图是时间序列与自身滞后的分布图,通常用来检验自相关性。如果序列中出现下图中的模式,那么说明该序列具有自相关性。如果没有出现这些模型,该序列很可能为随机白噪声。

    在下面太阳黑子区的例子中,随着滞后期的增长,图中的点越来越分散。

    from pandas.plotting import lag_plot
    plt.rcParams.update({'ytick.left' : False, 'axes.titlepad':10})
    
    # Import
    ss = pd.read_csv('https://raw.githubusercontent.com/selva86/datasets/master/sunspotarea.csv')
    a10 = pd.read_csv('https://raw.githubusercontent.com/selva86/datasets/master/a10.csv')
    
    # Plot
    fig, axes = plt.subplots(1, 4, figsize=(10,3), sharex=True, sharey=True, dpi=100)
    for i, ax in enumerate(axes.flatten()[:4]):
        lag_plot(ss.value, lag=i+1, ax=ax, c='firebrick')
        ax.set_title('Lag ' + str(i+1))
    
    fig.suptitle('Lag Plots of Sun Spots Area \n(Points get wide and scattered with increasing lag -> lesser correlation)\n', y=1.15)    
    
    fig, axes = plt.subplots(1, 4, figsize=(10,3), sharex=True, sharey=True, dpi=100)
    for i, ax in enumerate(axes.flatten()[:4]):
        lag_plot(a10.value, lag=i+1, ax=ax, c='firebrick')
        ax.set_title('Lag ' + str(i+1))
    
    fig.suptitle('Lag Plots of Drug Sales', y=1.05)    
    plt.show()
    

    在这里插入图片描述


    9 如何评估时间序列的可预测性?

    一个时间序列的模式越有规律,就越容易预测。可以用近似熵来量化时间序列的规律性和波动的不可预测性。近似熵越高,意味着预测难度越大

    另一个选择是样本熵。样本熵与近似熵类似,但在不同的复杂度上更具有一致性,即使是小型时间序列。例如,相比于“有规律”的时间序列,一个数据点较少的随机时间序列的近似熵较低,但一个较长的随机时间序列却有较高的近似熵。因此,样本熵更适于解决该问题。


    10 为什么要对时间序列做平滑处理?如果操作?

    对时间序列做平滑处理有以下几个用处:

    • 减少噪声影响,从而得到过滤掉噪声的、更加真实的序列。
    • 平滑处理后的序列可用作特征,更好地解释序列本身。
    • 可以更好地观察序列本身的趋势

    那么如果进行平滑处理呢?现讨论以下几种方法:

    • 取移动平均线
    • 做 LOESS 平滑(局部回归)
    • 做 LOWESS 平滑(局部加权回归)

    移动平均是指对一个滚动的窗口计算其平均值,该窗口的宽度固定不变。但你必须谨慎选择窗口宽度,因为窗口过宽会导致序列平滑过度。例如,如果窗口宽度等于季节长度,就会消除掉季节因素的作用。

    LOESS 是 Localized regression 的缩写,该方法会在每个点的局部近邻点做多次回归拟合。此处可以使用 statsmodels 包,你可以使用参数 frac 控制平滑程度,即决定周围多少个点参与回归模型的拟合。


    11 如何用格兰杰因果关系检验判断一个序列能否预测另一个?

    该检验基于一个假设,即 X 导致了 Y 的发生,那么基于 Y 的先前值与 X 的先前值得到的 Y 的预测值,要优于仅根据 Y 本身得到的预测值

    statsmodels 包提供了便捷的检验方法。它可以接受一个二维数组,其中第一列为值,第二列为预测因子。

    零假设为:第二列的序列与第一列不存在格兰杰因果关系。如果 P 值小于显著性水平 0.05,就可以拒绝零假设,从而知道 X 的滞后是起作用的。

    第二个参数 maxlag 表示该检验中涉及了 Y 的多少个滞后期。

    from statsmodels.tsa.stattools import grangercausalitytests
    df = pd.read_csv('https://raw.githubusercontent.com/selva86/datasets/master/a10.csv', parse_dates=['date'])
    df['month'] = df.date.dt.month
    grangercausalitytests(df[['value', 'month']], maxlag=2)
    

    参考资料:
    https://www.machinelearningplus.com/time-series/time-series-analysis-python/

    展开全文
  • Matlab时间序列分析

    万次阅读 多人点赞 2018-11-13 18:53:46
    在引入时间序列前,先介绍几个matlab函数 matlab中的gallery函数简析 Matlab 中的 gallery 函数是一个测试矩阵生成函数。当我们需要对某些算法进行测试的时候,可以利用gallery函数来生成各种性质的测试矩阵。其用法...

    时间序列分析需要解决的问题

    假设我们有这么一段数据(采样自移动公司某段时间的开户客户数目)
    下载地址,包含本实验用到的数据/资料以及全部代码,实验报告,引用请告知,关注并留言获取提取码,或者到我的csdn下载进行下载
    在这里插入图片描述
    我们可以发现,在数据中似乎有种明显的趋势规律,或许还存在一定的周期性。但是,我们却不好描述这种趋势或者说是周期,因为数据中存在一定的随机因素模糊了我们的判断。在进行时间序列数据分析时,我们大部分的目的都是想从中找出可以表达的规律以便于我们对未来做预测。
    在这里,我们把时间序列数据分解为以下三项:

    真实数据=趋势项+周期项+随机项
    趋势性:例如随时间变化的一次函数/多次函数/幂函数趋势等
    周期项:周期规律
    随机项:在这里可以理解为噪声,但是这些噪声也是存在规律的。
    

    不难发现,趋势性和周期项都不难描述,难以刻画的是随机项变化。而如果我们真正能预测随机项的变化,那么我们就基本上可以预测整个数据下一步的走势了。

    时间序列分析的步骤

    • 观测数据(均值,周期等)
    • 数据预处理
    • 平稳化–去趋势与去周期,剩余随机项
    • 用类似于回归/滑动平均的思想来 拟合随机项
      1.判断去趋势去周期后的数据是否平稳
      2.计算数据的自相关函数和偏相关函数
      3.根据自相关/偏相关函数性质决定选用什么模型来拟合随机项

    在这里插入图片描述
    4.模型定阶和拟合参数的求解
    5.模型检验

    • 模型最终实现功能:将周期项/趋势性/随机项结合,是一个关于时间变化的函数,这个函数可以在一定程度上对未来数据进行预测。

    如何实现每个步骤

    去趋势/去周期

    运用差分算子可以去趋势和周期,其中N阶差分用于具有N阶多项式趋势的数据去趋势,N步差分用于具有N步周期的数据去周期。

    偏相关/自相关函数的计算

    可以运用matlab自带函数,也可以自己编写。

    模型定阶

    可以参考AIC/BIC等定阶方法

    模型检验

    通过计算LB统计量和卡方统计量,判断拟合的误差项是否为白噪声序列,如果是,则检验通过。

    一个具体的案例 --移动开户数分析

    我们有一份移动公司两年间共计700余天的每日开户树,我们想探寻其中的随机因素。

    数据观察

    变化趋势

    在这里插入图片描述

       图一.原开户数据趋势变化的观察
    

    可以发现,开户数在局部区间的极值变化具有一定的周期性。但是数据太多且有点杂乱。
    为了更好的观察数据的规律,拟将数据进行滑动平均处理,在一定程度上更能看出规律性。并且,依据两年的月数天数变化,设置标记每月天数变化的向量,计算每月,每个季度,每年的平均开户数进行观察。计算并绘制得到下图:
    在这里插入图片描述

    图二.开户数据趋势变化周期的观察
    

    上图红色标记为2012年,蓝色标记为2013年。总的看来均有上升的趋势。图二左上角绘制的是2012和2013的开户数滑动平均观察(使用matlab smooth函数,滑动窗口设置为30),容易发现在年这个层面上趋势具有较大的相关性,一年显然可以是一个大周期的刻画。但是这个周期太大了,如果选取这个周期做去季节性,显然会损失大量数据。因此,我们想找到一个更小的周期来刻画,其实我们还可以发现图一中有明显的小起伏,经验证起伏平均间隔为30即一个月,那么一个月是否为一个周期呢?右上图是按月进行平均的统计值绘制,起伏貌似没有什么规律。但是还是可以作为参考。左下角图绘制的是季度的平均开户数变化,可以发现两个季度为一个周期,因为起伏在两个周期中进行。右下角图绘制的是年平均变化图,如我们所料,随着移动互联网的普及,2013年开户平均数比2012年高。

    综上观察,我们可以得出以下结论:
    1.一个月应该是一个最小周期,两个季度是一个中周期,一年为一个确定的大周期。
    2.开户人数有明显的上升趋势。
    3.由观测值我们可以认为存在异常值。

    数据预处理

    2.1缺失值处理
    数据中不存在缺失值
    2.2异常值处理
    拟采用三倍标准差即拉以达法则筛选异常值,考虑到时间序列的短期影响性,用其周围值平均代替。

    去趋势与去周期

    进行了简单的观察和数据预处理后,我们开始正式进行时间序列平稳化的操作。考虑到有两年,我们从单独考虑2012,2013,然后集中考虑2012和2013来进行处理和检验。
    3.1去趋势
    由于数据变化有一定的集中性,且有滑动平均看来近似可以用一次函数拟合。所以采用一次差分的方法去除数据的趋势项。这里采用diff函数实现一次差分。去趋势后,我们进行了如下的观察:

    在这里插入图片描述

    在这里插入图片描述

    图三.周期为一个月的检验图
    

    我们发现相邻极值点间的时间差距十分接近于一个月,这让我们欣喜万分。

    3.2去周期
    由以上的观察数据和去趋势后数据综合考虑,采用一个月作为周期是个不错的选择。这里自定义d步差分函数对去趋势后的数据进行去周期的处理。
    经过3.1和3.2,我们得到如下结果:
    在这里插入图片描述在这里插入图片描述

    图四.2012年去趋势和去周期比较图 图五.2013年去趋势和去周期比较图

    在这里插入图片描述

    图六.2012-2013年去趋势和去周期比较图
    

    肉眼看来,效果还算不错,从均值线看来均值均十分接近于0。为了更客观的看待平稳的效果,采用自相关图还进行平稳性的验证。

    自相关性平稳性检验结果如下:
    在这里插入图片描述

    图七.自相关图检验平稳性
    

    从自相关图看来,经过去趋势和去周期后的数据的确是平稳的。 经过以上步骤,我们就成功的将非平稳的时间序列转换为了平稳的时间序列,得到了随机误差项的表示。

    截至目前以上代码总结:

    %% 计算月平均,季度平均,年平均意图找寻周期规律
    % 输入开户数据
    % 输出各有用的统计量
    function [months_mean,seasons_mean,years_mean,month_starts,month_ends,season_starts,season_ends] = Calu_mean(open_nums)
    t = 1:length(open_nums);
    nums = open_nums;
    month_days = [0,31,29,31,30,31,30,31,31,30,31,30,31,31,28,31,30,31,30,31,31,30,31,30,31]; 
    season_days = [0,sum(month_days(2:4)),sum(month_days(5:7)),sum(month_days(8:10)),sum(month_days(11:13))...
        ,sum(month_days(14:16)),sum(month_days(17:19)),sum(month_days(20:22)),sum(month_days(23:25))];
    % calu the month_mean nums
    months_mean = zeros(1,24);
    month_starts = zeros(1,24);
    month_ends = zeros(1,24);
    for ii = 1:length(month_days)-1
        start_day = sum(month_days(1:ii))+1;
        end_day = sum(month_days(1:ii+1));
        month_starts(ii) = start_day;
        month_ends(ii) = end_day;
        months_mean(ii) = mean(nums(start_day:end_day));
    end
    
    %% 去趋势和周期函数
    % 输入预处理后的时序数据向量vector和观察数据变化人为认定的周期T
    % 输出去趋势后的向量数据detrend_data和去趋势去周期后的向量数据detrend_deT_data
    function [detrend_data,detrend_deT_data] = Detrend_plot(vector,T)
        subplot(3,1,1)
        plot(1:length(vector),vector,'k')
        hold on
        plot(0:length(vector),mean(vector)*ones(1,length(0:length(vector))),'r')
        text(1,mean(vector),'均值')
        legend('趋势','均值')
        title('原开户数据趋势')
        % 去趋势项
        detrend_data = diff(vector,1);
        subplot(3,1,2)
        plot(1:length(detrend_data),detrend_data,'k')
        hold on
        plot(0:length(detrend_data),mean(detrend_data)*ones(1,length(0:length(detrend_data))),'r')
        legend('趋势','均值')
        title('去除趋势项后趋势')
        % 去季节项
        detrend_deT_data = zeros(length(detrend_data),1);
        for i=1:length(detrend_data)
            if(i<=T)
                detrend_deT_data(i) = [];
            else
                detrend_deT_data(i) = detrend_data(i)-detrend_data(i-T);
            end
        end
        subplot(3,1,3)
        plot(1:length(detrend_deT_data),detrend_deT_data,'k')
        hold on
        plot(0:length(detrend_deT_data),mean(detrend_deT_data)*ones(1,length(0:length(detrend_deT_data))),'r')
        legend('趋势','均值')
        title('去除趋势项后和季节项后趋势')
    end
    
    %% 主函数脚本,运行该脚本将输出以上所有结果
    clearvars
    %% 读取数据并预处理
    open_nums = xlsread('移动通知户开户数.xlsx',1,'B2:B732');
    %% 计算月平均,季度平均,年平均并绘图找周期规律,并处理异常值
    figure(1),[months_mean,seasons_mean,years_mean,month_starts,month_ends,season_starts,season_ends] = Calu_mean(open_nums);
    [m,n] = find(abs(open_nums-mean(open_nums))>2*std(open_nums));
    open_nums(m) = mean(open_nums);
    %% 尝试进行去趋势和去周期,考虑到存在多重周期,将年份分开按照月为周期去除周期(其中每月为多少天根据两年分别不同)
    data_2012 = open_nums(1:366);
    data_2013 = open_nums(367:end);
    figure(2),[detrend_data2012,detrend_deT_data2012] = Detrend_plot(data_2012,30);% 2012开户数变化趋势
    figure(3),[detrend_data2013,detrend_deT_data2013] = Detrend_plot(data_2013,31);% 2013开户数变化趋势
    figure(4),[detrend_data,detrend_deT_data] = Detrend_plot(open_nums,30);% 2012-2013开户数变化趋势
    figure(5),subplot(411),autocorr(open_nums),title('原数据自相关图')% 自相关图分析
    subplot(412),autocorr(detrend_data),title('2012年随机项自相关图')
    subplot(413),autocorr(detrend_deT_data2013),title('2013年随机项自相关图')
    subplot(414),autocorr(detrend_deT_data2013),title('2012-2013年随机项自相关图')
    clearvars month_starts month_ends season_starts season_ends data_2012 data2013
    

    截至目前,我们编写的代码所完成的工作有:
    1.观察数据
    2.数据预处理
    3.去趋势,去周期
    4.计算了偏相关函数和自相关函数

    我们还需要进行的工作是:
    1.判断处理后的随机项是否平稳,判断自相关/偏相关函数的拖尾性/截尾性,以确定采用什么模型来拟合随机项。去均值使该时间序列变为零均值平稳过程。
    在这里插入图片描述
    2.确定模型后,确定采用几阶进行拟合,并求解参数
    3.模型检验

    ******* 休息好了吧,我们继续

    判断去趋势和去周期后数据的平稳性

    在这里插入图片描述在这里插入图片描述

     图一  处理后的数据               图二  处理后数据的20期自相关系数
    

    自相关系数定义为:
    在这里插入图片描述
    由图二观察知道,当|h|>1时,自相关系数为近似为0.认为自协方差函数(自相关函数)有“截尾性”.可以使用MA(1)模型。为了考虑能否使用ARMA模型,进一步进行讨论。

    转换为零均值平稳序列

    采用mean函数计算该序列的均值,处理如下:
    随机变量总体X 的样本的均值估计量:
    在这里插入图片描述
    在这里插入图片描述

    转换前后序列变化如下:
    在这里插入图片描述在这里插入图片描述

      图三 非零均值平稳序列                图四 零均值平稳序列 
    

    3.无偏/有偏自相关函数的估计
    自协方差函数(或自相关函数)的样本估计量通常有两种类型:
    在这里插入图片描述(3.3.1)
    在这里插入图片描述 (3.3.2)
    关于两种估计量,前人曾经做过比较,G.M.Jenkins 结论:
    1.估计量( 3.3.1) 更满意。
    2.其中(3.3.1)是自相关函数的有偏估计,(3.3.2)是自相关函数的无偏估计。
    3.(3.3.1)式定义的估计量收敛于零的速度更快。
    4.由(3.3.1)式定义的样本自协方差函数矩阵(或相关矩阵)是正定的

    根据两种自相关系数的计算公式,我们可以自己编写函数计算自相关函数如下
    其中函数输入为零均值的平稳序列,输出为一组无偏估计的自相关系数,一组为有偏的自相关系数。完成自己的函数编写后,选择计算前面30期自相关系数与matlab自带autocorr进行比较。
    自定义计算无偏和有偏两种自相关系数代码:

    function [unbiased_autocorr,biased_autocorr] = my_autocorr(time_vector)
    %% 函数定义,计算无偏和有偏两种自相关系数
    n = length(time_vector);
    unbiased_autocorr = zeros(n,1);
    biased_autocorr = zeros(n,1);
    % 延迟期数从0到n
    for h = 0:n
        omega = 0;
        % 累加从t=1到n-h
        for t = 1:n-h
            omega = omega+time_vector(t)*time_vector(t+h);
        end
        unbiased_autocorr(h+1) = omega/(n-h);
        biased_autocorr(h+1) = omega/n;
    end
    % 自相关函数转换为自相关系数
    unbiased_autocorr = unbiased_autocorr/unbiased_autocorr(1);
    biased_autocorr = biased_autocorr/biased_autocorr(1);
    %% 绘图比较
    n_compare = 30;
    figure()
    subplot(311) 
    plot(0:n_compare-1,unbiased_autocorr(1:n_compare),'ro','MarkerFaceColor','r','markersize',4)
    hold on,plot(0:n_compare-1,zeros(1,n_compare),'b')
    title('无偏自相关系数变化'),xlabel('延迟天数'),ylabel('自相关系数')
    subplot(312) 
    plot(0:n_compare-1,biased_autocorr(1:n_compare),'ro','MarkerFaceColor','r','markersize',4)
    hold on,plot(0:n_compare-1,zeros(1,n_compare),'b')
    title('有偏自相关系数变化'),xlabel('延迟天数'),ylabel('自相关系数')
    subplot(313) 
    autocorr(time_vector,n_compare),title('内置autocorr结果')
    

    结果如图:
    在这里插入图片描述

         图五  前30期自相关系数比较图
    

    由图,计算结果与内置函数结果十分接近,验证了自己构建的函数的准确性。
    在这里插入图片描述

    图六 无偏/有偏自相关系数的收敛速度比较
    

    可以看出。该图验证了第三个结论,即有偏估计定义的估计量收敛于零的速度更快,且从收敛趋势看来效果更好。

    4.偏相关系数的估计和比较
    计算偏相关系数有两种方法,直接求解和递推求解。这里考虑递推求解。运用3求解得到的自相关系数一步步递推得到所有的偏相关系数。
    4.1 直接求解
    在这里插入图片描述(3.4.1)
    由最优线性估计的计算公式,应用矩阵按行(列)展开式的性质可得
    在这里插入图片描述(3.4.2)
    4.2 递推求解
    在这里插入图片描述(3.4.3)
    递推过程:
    在这里插入图片描述

           图七 递推求解过程
    

    理解以上递推思想后,利用3计算出来的自相关系数,编写函数如下:
    其中输入为自相关系数和之前的零均值序列。前者用于递推,后者用于利用matlab自带函数parcorr对自定义的函数计算值进行检验。
    代码:

    function [customed_parcorr] = my_parcorr(time_vector,auto_corr)
    %% 函数定义
    n = length(auto_corr);
    customed_parcorr = zeros(n,n);
    customed_parcorr(1,1) = auto_corr(2);
    for k = 1:n-2
        sum1 = 0;
        sum2 = 0;
        for j = 1:k
            sum1 = sum1 + auto_corr(k+2-j)*customed_parcorr(k,j);
            sum2 = sum2 + auto_corr(j+1)*customed_parcorr(k,j);
        end
        customed_parcorr(k+1,k+1) = (auto_corr(k+2)-sum1)/(1-sum2);
        for j = 1:k
            customed_parcorr(k+1,j) = customed_parcorr(k,j) - customed_parcorr(k+1,k+1)*customed_parcorr(k,k+1-j);
        end
    end
    %% 绘图比较
    a = diag(customed_parcorr);
    subplot(211)
    compare_n = 40;
    parcorr(time_vector,compare_n)
    subplot(212)
    plot(1:compare_n,a(1:compare_n),'ro','markerfacecolor','r','markersize',4)
    hold on,plot(1:compare_n,zeros(1,compare_n),'b')
    ylim([-0.5,1])
    title('自定义偏相关计算值')
    xlabel('延迟天数'),ylabel('偏相关系数')
    

    结果如图八所示:
    在这里插入图片描述

    图八 自定义偏相关系数计算值与parcorr前30期计算值比较
    

    由比较图看来,计算结果相当。
    由图八可知,偏相关系数具有拖尾性。

    截至目前,除了上面给过的代码外,新的部分计算代码给出如下:

    %% 主函数运行脚本
    clearvars
    %% 读取数据并预处理
    open_nums = xlsread('移动通知户开户数.xlsx',1,'B2:B732');
    %% 计算月平均,季度平均,年平均并绘图找周期规律,并处理异常值
    figure(1),[months_mean,seasons_mean,years_mean,month_starts,month_ends,season_starts,season_ends] = Calu_mean(open_nums);
    [m,n] = find(abs(open_nums-mean(open_nums))>2*std(open_nums));
    open_nums(m) = mean(open_nums);
    %% 尝试进行去趋势和去周期,考虑到存在多重周期,将年份分开按照月为周期去除周期(其中每月为多少天根据两年分别不同)
    data_2012 = open_nums(1:366);
    data_2013 = open_nums(367:end);
    figure(2),[detrend_data2012,detrend_deT_data2012] = Detrend_plot(data_2012,30);% 2012开户数变化趋势
    figure(3),[detrend_data2013,detrend_deT_data2013] = Detrend_plot(data_2013,31);% 2013开户数变化趋势
    figure(4),[detrend_data,detrend_deT_data] = Detrend_plot(open_nums,30);% 2012-2013开户数变化趋势
    figure(5),subplot(411),autocorr(open_nums),title('原数据自相关图')% 自相关图分析
    subplot(412),autocorr(detrend_data),title('2012年随机项自相关图')
    subplot(413),autocorr(detrend_deT_data2013),title('2013年随机项自相关图')
    subplot(414),autocorr(detrend_deT_data),title('2012-2013年随机项自相关图')
    clearvars month_starts month_ends season_starts season_ends data_2012 data2013
    %% 平稳序列去均值变为零均值平稳序列
    zeromean_detrend_deT_data = detrend_deT_data-mean(detrend_deT_data);
    [unbiased_autocorr,biased_autocorr] = my_autocorr(zeromean_detrend_deT_data);
    %% 计算偏相关函数
    [customed_parcorr] = my_parcorr(zeromean_detrend_deT_data,biased_autocorr);
    

    ************** 休息一会,我们继续

    处理后数据适用的模型判别

    观察平稳序列的自相关函数和偏相关函数,通过是否拖尾/结尾判断使用的模型。
    在这里插入图片描述

    图*1 模型判别流程图
    

    通过上面的分析,自相关系数和自偏相关系数均呈现拖尾性,故使用ARMA模型进行进一步分析。

    3.ARMA模型的参数估计和定阶
    拟函数法进行求解,求解步骤如下:

    3.1初定阶数p和q,运用Yule-Walker求解逆函数I
    3.2求解MA滑动平均系数
    3.3求解AR自回归系数
    3.4在求出所有参数后,利用极大似然法残差平方和求方差,算出AIC/BIC指标
    3.5重复不走3.1-3.4,在一定大范围内找出极小BIC值对应的阶数p,q
    

    4.参数检验和定阶的进一步修正
    通过模型的显著性检验和参数的显著性检验,在参数估计和定阶最小BIC情况下对应的结束p,q附近进行,最终找到最佳的ARMA阶数为ARMA(1,1),参数求解为
    在这里插入图片描述
    经检验残差序列为白噪声,认为模型效果可以接受。
    1.适用模型识别
    首先,我们需要做的是根据动态数据,从各类模型族中选择与实际过程相吻合的模型;
    即处理后的平稳序列应该用什么模型来表示最为合理。经过前期的预处理,去趋势和去周期,得到平稳时间序列的自相关系数,偏相关系数如下图1.1。

    在这里插入图片描述在这里插入图片描述

                    图1.1 平稳序列的自相关/自协相关系数
    

    不难发现,该平稳序列的自相关系数和自偏相关系数均有拖尾的性质,且自相关拖尾阶数较高分布在32-35阶左右,而自偏相关系数拖尾阶数较低,分布在12-20之间。综上,拟采用ARMA模型对其进行拟合。
    2.ARMA模型的参数估计和定阶
    确定好适用的模型之后,下一步我们要进行的是对确定的模型类别,定出模型中的阶数p 和 q ,这里采用先大范围搜索,在搜索的结束中求的模型的参数估计——自回归系数和滑动平均系数,并由此求出拟合残差的方差。最后,计算每个阶数下对应的AIC/BIC进行比较,这里适用最小BIC值对应的阶数和参数估计做下一步的分析。

    2.1 方法的选取
    课件4.3关于ARMA模型的参数估计中介绍了两种估计参数的方法,对这两种方法的理解如下:
    *ARMA模型参数的矩估计:通过Y-W方程求解出自回归系数,利用MA(q)模型参数的线性迭代法或牛顿—拉裴森算法求得滑动平均系数,最后合并可得拟合模型。这种方法 矩估计方法精度较低,通常只能作为其他迭代参数估计算法的初值。
    *模型参数的初估计法—逆函数法:利用可逆域的逆转形式进行求解。这种方法仅涉及到线性方程的求解,计算简单便于实现。
    综合上述两种方法,考虑到逆函数法对于低阶或高阶模型都很有效且求解方便。这里采用逆函数法对模型参数进行估计,逆函数法求解步骤如下
    逆函数法求解步骤

    1.令P(*) = max(p,q)+q,用y—w方法估计AR(p*)的 p* 个自回归参数,作为
    相应的ARMA( p, q)模型的前 p*个逆函数值.
    2.求解线性方程组得到滑动平均系数的估计。
    3.将前两步的结果代入课件4.3.6公式,求得自回归系数的估计。
    

    模型参数和求解和定阶

    2.2.1模型的定阶指标
    *最终预报误差准则(PDE)
    在这里插入图片描述
    其中,方差依赖于k, 其大小反映了模型与数据的拟合程度。第一个因子随 k 增大而增大, 放大了残差方差的不确定性影响.
    *最小信息准则(AIC)
    在这里插入图片描述
    AIC准则可应用于ARMA模型及其他统计模型(如多项式回归定阶).能在模型参数极大似然估计的基础上,对于ARMA(p, q)的阶数和相应的参数,同时给出一种最佳估计.
    *BIC准则
    在这里插入图片描述
    AIC准则避免了统计检验中由于选取置信度而产生的人为性,为模型定阶带来许多方便. 但AIC方法未给出相容估计.
    **注意:AIC阶数估计一般偏高,BIC阶数估计偏低。这里适用BIC准则进行初始阶数的定阶,取BIC值最小下对应的阶数p,q即为初始定阶。
    2.2.2 模型的求解
    观察之前的自相关系数图和自偏相关系数图,初定ARMA模型初始阶数的选取p和q均在1-20阶内,利用逆函数求解并算的每个p,q下参数对应的BIC值,取最小BIC对应的p,q即为初定的阶数。不同怕p,q下对应的BIC值如下图。
    在这里插入图片描述

    图2.2.2 BIC值的变化
    

    由BIC值的变化我们可以看出,在取ARMA(1,1)时BIC值最小。此时ARMA(1,1)模型建立如下:
    在这里插入图片描述

    定阶的进一步修正-模型的检验

    确定好最小BIC准则下的参数估计值和模型阶数p和q后,我们知道——以上各个问题是相互关联的,需整体进行系统化的模型优化. 而根据已有的一组样本数据建立模型,对模型阶数和参数作出判断和估计,是重要的两部分工作.所以最后,我们需要在已经初定的阶数周围小范围内进行搜索和检验,通过模型的有效性检验的参数的显著性检验进一步优化模型,寻找最优的阶数和参数估计值。统计模型只是对生成观测数据真实过程的近似,在模型拟合后,还需要进行模型的有效性检验,及检验拟合模型对序列中信息的提取是否充分。模型的有效性检验即对残差的白噪声检验。
    观察残差序列及其自相关系数如下图所示:

    在这里插入图片描述

    图3.1 残差序列自相关系数
    

    可以发现序列间没有相关性,不具有拖尾和截尾的性质,可以初步判断改模型是可靠的。为了进一步验证残差序列是否为白噪声序列,下面进行LB统计量的检验结果如下:
    延迟期数 LB统计量 卡方统计量

    在这里插入图片描述
    由LB统计量小于对应的卡方统计量,认为该序列可以认为是白噪声序列。模型效果较好,模型检验通过。
    到此,我们所有工作已经完毕。如果有什么疑问,欢迎留言,我会一一解答。

    新部分的代码:
    1主函数 main1.m
    实现功能:读取处理后的平稳数据进行建模,通过观察自相关/自偏相关系数确定选用ARMA模型,在一定确定的阶数范围内求模型的参数估计,根据最小BIC准则初定p,q的选择,最后通过模型的有效性检验和参数的显著性检验确定最终的参数估计和阶数。

    % 读取处理后的平稳时间序列
    
        processed_data = csvread('detrend_deT_data.csv');
        figure(1)
        autocorr(processed_data,60)
        figure(2)
        parcorr(processed_data,25)
        data_autocorr = autocorr(processed_data,600);
        % 自相关系数向量转换为自相关系数矩阵
        autocorr_matrix = vec2mat(data_autocorr);
        max_p = 15;
        max_q = 20;
        AICs = ones(max_p,max_q);
        BICs = ones(max_p,max_q);
        N = 600;
        for p=1:20
            for q = 1:20
                % 计算I
                P = max(p+q)+q;
                I = autocorr_matrix(1:P,1:P)\autocorr_matrix(1,2:P+1)';
                % 计算MA部分的滑动平均系数
                seta = calu_MA_parameter(I,p,q);
                % 计算AR部分的自回归系数
                fine = calu_AR_parameter(seta,I,p,q);
                % 计算方差
                residual_var = calu_var(p,q,fine,seta,N,processed_data);
                % 计算AIC和BIC并存储
                [AICs(p,q),BICs(p,q)] = calu_AIC_BIC(p,q,N,residual_var);
            end
        end
    

    % 找出最小BIC下的阶数p,q
    [min_BIC_p,min_BIC_q] = find(BICs == min(BICs));
    2.MA滑动平均系数的求解函数
    [seta]=calu_MA_parameter(I,p,q)
    函数输入:逆函数{ Ij }的估计I,给定的自回归和滑动平均阶数p和q
    函数输出:滑动平均系数的估计值seta

    function [seta]= calu_MA_parameter(I,p,q)
    P = max(p,q);
    I_matrix = zeros(q,q);
    for ii = 1:q
        I_matrix(ii,:) = I(P+ii-1:-1:P+ii-q);
    end
    seta = I_matrix\I(P+1:P+q);
    end
    

    3.AR自回归系数的求解函数
    [fine] = calu_AR_parameter(seta,I,p,q)
    函数输入:逆函数{ Ij }的估计I,给定的自回归和滑动平均阶数p和q以及滑动平均系数seta
    函数输出:自回归系数的估计值fine

    function [fine] = calu_AR_parameter(seta,I,p,q)
    fine = ones(1,p);
    for ii = 1:p
        if ii>q
            seta(ii) = 0;
        end
        fine(ii) = seta(ii)+I(ii);
        for jj = 1:ii-1
            if (ii-jj)>0
                fine(ii) = fine(ii)-seta(jj)*I(ii-jj);
            end
        end
    end
    

    4.残差方差的计算函数
    [residual_var]= calu_var(p,q,fine,seta,N,processed_data)
    函数输入:逆函数{ Ij }的估计I,给定的自回归和滑动平均阶数p和q以及滑动平均系数seta,自回归系数fine,处理后的平稳序列。
    函数输出:残差的方差

    function [residual_var] = calu_var(p,q,fine,seta,N,processed_data)
    error = zeros(1,N);
    for ii = max(p,q)+1:N
        error(ii) = [1,-fine]*processed_data(ii:-1:ii-p)+seta'*error(ii-1:-1:ii-q)';
    end
    residual_var = (error*error')/N;
    end
    

    5.AIC/BIC指标的计算函数
    [AIC,BIC] = calu_AIC_BIC(p,q,N,residual_var)
    函数输入:给定的自回归和滑动平均阶数p和q,样本个数N以及残差方差。
    函数输出:不同阶数下对应的AIC/BIC指标。

    function [AIC,BIC] = calu_AIC_BIC(p,q,N,residual_var)
    AIC = log(residual_var)+2*(p+q+1)/N;
    BIC = log(residual_var)+log(N)*(p+q+1)/N;
    End
    

    6.计算LB统计量函数

    lb = calu_LB(data, m)
    函数输入:残差序列,延迟期数
    函数输出:LB卡方统计量
    function lb = calu_LB(data, m)
    N = 25;
    power_autocorr = 0;
    data_autocorr = autocorr(data,25);
    for k=1:m
        power_autocorr = power_autocorr+data_autocorr(k+1)^2/(N-k);
    end
    lb = N*(N+2)*power_autocorr;
    End
    

    7.自相关系数向量转为矩阵 vec2mat
    [autocorr_matrix] = vec2mat(autocorr_vector)
    函数输入:自相关系数函数输出:LB卡方统计量

    function [autocorr_matrix] = vec2mat(autocorr_vector)
    n = length(autocorr_vector);
    autocorr_matrix = zeros(n);
    for ii = 1:n
        autocorr_matrix(ii,ii:n) = autocorr_vector(1:n-ii+1);
    end
    autocorr_matrix = triu(autocorr_matrix)+tril(autocorr_matrix',-1);
    end
    

    yqol 编写不易,欢迎关注。

    展开全文
  • 时间序列预测分析方法(一):相关分析

    万次阅读 多人点赞 2019-06-18 14:35:11
    针对特定的预测问题,只是拥有数据还不够,想要从纷繁复杂的数据关系中挖掘出可用于预测的规律或模式,还得运用恰当的分析方法。比如聚类分析,恰当地选择聚类算法,可以按维度将数据适当地分群,根据各类的特征制订...

            针对特定的预测问题,只是拥有数据还不够,想要从纷繁复杂的数据关系中挖掘出可用于预测的规律或模式,还得运用恰当的分析方法。比如聚类分析,恰当地选择聚类算法,可以按维度将数据适当地分群,根据各类的特征制订营销计划或决策,抑或是根据各类不冋规律建立起更有针对性的预测模型;还有常用的关联分析,可以从事物的历史数据中挖掘出变化规律有指导性地对未来进行预测,如此等等。本内容将分别介绍常用的分析方法,本节介绍相关分析。

    目录

    1. 自相关分析            2. 偏相关分析           3. 简单相关分析           4. 互相关分析           5. 典型相关分析


           相关关系是一种与函数关系相区别的非确定性关系,而相关分析就是研究事物或现象之间是否存在这种非确定性关系的统计方法。相关分析按处理问题的不同,通常可分为自相关分析、偏相关分析、简单相关分析、互相关分析及典型相关分析。其中:

    • 自相关分析、偏相关分析:适用于分析变量自身的规律;
    • 简单相关分析:通常可分析任意两个等长数列间的相关性;
    • 互相关分析:允许在一定的间隔下进行简单相关分析;
    • 典型相关分析:适用于分析两组变量的相关性。

    变量与变量的关系通常有两种:

    • 函数关系:表示确定的非随机变量之间的关系,可以用表达式具体地定义。比如,位移与速度就是确定的函数关系。
    • 非确定性关系:表示非随机变量与随机变量的关系,也就是说当给定一个变量时,另一个变量是随机的,它不能由具体的表达式来定义。比如,身高和体重就无法用一个确定的函数关系来表达,但是根据经验,身高比较高的人,其体重也比较重,说明两者是相关关系

    1. 自相关分析

            自相关是指同一时间序列在不同时刻取值的相关程度,假设有时间序列X_{t}, t=1,2,3,...,则在时刻t和滞后n阶t+n之间的相关即为n阶自相关,其定义如下: 

                                

    其中,函数f为计算相关系数的函数,可通过上式计算滞后n阶自相关系数的值。这里使用 airmiles时序数据来分析时间序列的自相关性,该数据集记录的是从1937年到1960年美国商业航空公司的飞机里程数据,如图1-1所示。 

                           

                                                      图 1-1 从1937年到1960年美国商业航空公司的飞机里程趋势

    在R语言中,可用acf函数分析时间序列的自相关性,该函数定义及参数说明如下表: 

    acf(airmiles,type='correlation',lag.max=10)
    
    # airmiles:美国1937-1960年客运里程营收(实际售出机位乘以飞行哩数)

    效果图如下: 

                       

    如图可见,滞后阶数为0时,相关系数为1,随着滞后阶数的增加,相关系数逐渐减弱,并趋于稳定。 

    2. 偏相关分析 

           在分析当期与前n期指标的相关性时,1期到n-1期之间的指标会对相关性的分析有影响,于是将这些影响指标去掉后,再进行的相关分析叫作偏相关分析。

            由于自相关性分析的是时间序列X_{t}, t=1,2,3,...,在时刻tt+n之间取值的相关性程度,其值是未在限定X_{t+1}\sim X_{t+n-1}取值的情况下进行计算的,所得的自相关系数多少会受X_{t+1}\sim X_{t+n-1}取值的影响。为了更加真实地计算自相关关系数值,需要在限定其他值的情况下进行计算,这就是所谓的偏相关,其定义如下:

            

          其中pf函数是求解X_{t}X_{t+n}在排除X_{t+1}\sim X_{t+n-1}因素影响的情况下的偏相关系数。同时,对X_{s}X_{t}在限定X_{k}的情况下,其偏相关系数定义如下: 

                                           

         在R语言中,可直接使用pacf函数分析时间序列的偏相关性,该函数定义及参数说明如表3-1-2所示 :

             

    pacf(airmiles,lag.max=10)

     效果图如下: 

     

           如图所示,最小为1阶滞后,对应值为0.876,与对应的1阶自相关系数相等,随着滞后阶数的增加(大于2阶),偏相关系数一直较小并且稳定。

    1阶滞后:相当于X_{t}X_{t+1},因此其自相关系数和偏自相关系数相等。

    3. 简单相关分析

           相关关系是一种非确定的关系,就好像身高与体重的关系一样,它们之间不能用一个固定的函关系来表示。而相关分析就是研究这种随机变量间相关关系的统计方法。此处,主要探讨不同特征对研究对象的相关性影响。常见进行相关分析的方法,主要有散点图和相关图。

    3.1 散点图

           散点图就是数据点在直角坐标系上的分布图,通常分为散点图矩阵和三维散点图。其中散点矩阵是由变量两两组合由数据点分布图构成的矩阵,而三维散点图就是从所有变量中选择三个变量进行绘制,进一步在三维空间里观察数据的形态。

    (1)散点图矩阵

           常用于绘制散点矩阵的方法是R语言中:

    • car包中的 scatterplotMatrix函数;
    • graphics包中的pais函数。

           这里以R语言自带的iris数据集为例,说明分析鸢尾花的四个指标的相关关系,绘制散点图矩阵,代码如下。

    # 图3-1
    pairs(~Sepal.Length+Sepal.Width+Petal.Length+Petal.Width,data=iris, 
       main="Simple Scatterplot Matrix")
    
    # 图3-2
    library(car)
    scatterplotMatrix(~Sepal.Length+Sepal.Width+Petal.Length+Petal.Width|Species, 
           data=iris)
    
    # 图3-3
    library(rgl)
    scatter3d(iris$Sepal.Length, iris$Petal.Length, iris$Petal.Width)

                

                                                                                          图 3-1

           如图3-1所示为所有变量的两两组合的散点图矩阵,每个散点图中呈现的是任意两变量的数据点,可通过数据点的分布,了解变量之间的相关性。此图中 Petal.Length与 Petal.Width对应的散点图比较接近线性,说明这两个变量间的相关性较强。

               

                                                                                                     图 3-2

           如代码所示,通过“|”指定了分组的变量,这里使用鸢尾花的种类进行分组。如图3-2所示对角线上的图形表示各个变量在不同鸢尾花类型下的分布情况;其他图形分别用不同颜色为数据点着色,同时实线表示对数据点的拟合线,上下虚线表示浮动范围,直接表示线性拟合线。根据该图可以更进一步地知道不同类型鸢尾花各变量的相关关系,以及线性及非线性的变化规律

    (2)三维散点图

            常用于绘制三维散点图的方法是:

    • car包的 scatter3d函数

            通过该函数绘制的图形,可以进行转动,以便切换不同角度,查看数据在三维空间的分布情况。通过使用 scatter3d函数绘制鸢尾花的 Sepal.Length、 Petal.Length、 Petal.Width这三个指标在三维空间的散点图,如下所示,该函数为三维空间中的点拟合了线性平面,通过旋转可以更直观地观察数据的分布规律。

                                             

                                                                                                图 3-3

    3.2 相关图

            所谓相关图是基于变量间的相关关系大小,通过可视化方式反应不同变量组合间相关关系的差异的图形。可以把相关图分为相关矩阵图、相关层次图。

    (1)相关矩阵图

           R语言中,绘制相关矩阵图的包,主要有两个,分别为:

    • corrgram包中的 corrgram函数;

            panel:统一设置非对角的面板,

            lower.panel和 upper panel:分别设置相关矩阵下三角和上三角的面板形状,

            通常可以设置:panel.pie(饼图)、 panel.shade(阴影图)、 panel.ellipse(椭圆图)、panel.pts(散点图)

    library(corrgram)
    
    #1、设置排序处理
    corrgram(mtcars,order=TRUE)
    
    #2、设置上下三角面板形状
    corrgram(mtcars,order=TRUE,lower.panel=panel.shade,upper.panel=panel.pie)
    
    #3、只显示下三角部分
    corrgram(mtcars,order=TRUE,lower.panel=panel.shade,upper.panel=NULL)
    
    #4、调整面板颜色
    corrgram(mtcars,order=TRUE,lower.panel=panel.shade,upper.panel=panel.pie,
             col.regions=colorRampPalette(c("darkgoldenrod4","burlywood1","white",
             "darkkhaki","darkgreen")))

     

                                                                    图 3-4:以上程序1、2、3、4效果图

    • complot包中的 complot函数

             method:可以取7个值,分别为 circle(圆形,默认) 、square(方形) 、ellipse(椭圆)、number(数字)、pie(饼图) 、shade(阴影) 、color(热力图)。

            type:用来设置相关矩阵图的显示区域,通常取full全部,默认)upper(上三角)、lower(下三角)

            col:可通过 colorRampPalette函数向col参数赋值来设置颜色

            order:用来指定相关系数变量的排序方式,通常可取 original(原始顺序,默认),AOE(特征向量的角顺序,FPC(第一主成分的顺序), hclust(层次聚类的顺序)以及 alphabet(字母表顺序)

    library(corrplot)
    
    #1、设置method=color绘制热力矩阵图
    corrplot(cor(mtcars), method="color", order = "AOE",tl.col="black",tl.srt=45,
             addCoef.col="black",col=colorRampPalette(c("#7F0000","red","#FF7F00",
             "yellow","white", "cyan", "#007FFF", "blue","#00007F"))(20))
    
    #2、绘制上下三角及不同色彩的相关矩阵图
    library(RColorBrewer)
    par(mfrow=c(2,2))
    corrplot(cor(mtcars),type="lower")
    corrplot(cor(mtcars),type="lower",order="hclust", col=brewer.pal(n=8,name="RdYlBu"))
    corrplot(cor(mtcars),type="upper",order="AOE", col=c("black","white"),bg="lightblue")
    corrplot(cor(mtcars),type="upper",order="FPC", col=brewer.pal(n=8, name="PuOr"))
    par(mfrow=c(1,1))

        

                                                                         图 3-5:以上程序1、2效果图 

    (2)相关层次图

          层次图,通过计算样本间的距离来判断各样本是否属于同一类。通过将相关系数转化为距离度量,进行系统聚类,旨在分析各变量的相关关系及组合影响情况。通常有四种方法将相关系数转化为相异性度量d

     使用mtcars数据集,进行系统聚类,如下: 方式一:hclust() 

    • d=1-r
    • d=(1-r)/2
    • d=1-|r|
    • d=\sqrt{1-r^{2}}
    d<-sqrt(1-cor(mtcars)^2)
    hc<-hclust(as.dist(d))
    plot(hc)
    rect.hclust(hc,k=3)  # 添加聚类分类矩形,如分为3类

     效果如下图所示: 

                           

    方式二:pvclust(),通过p值来评估层次聚类的不确定性,表明了影响显著的聚类,显著的类越多,效果越好。

    参数:method.dist() 指定用于计算距离的方法,设定为correlation时,按 d=1-r 计算;设定为abscor时,按 d=1-|r| 计算。

    library(pvclust)
    cluster.bootstrap <- pvclust(mtcars, nboot=1000, method.dist="correlation")
    plot(cluster.bootstrap)
    pvrect(cluster.bootstrap)  # 添加聚类分类矩形

     效果如下图所示:

                         

    • AU值:Approxiamtely Unbiased,通过多尺度自助重抽样法计算而来;
    • BP值:Bootstrap Probability,通过普通自助重抽样法计算而来 

    4. 互相关分析

          与自相关不同,互相关是指两个时间序列在做任意两个不同时刻的相关程度,假设有时间序列X_{t}, t=1,2,3,...Y_{t}, t=1,2,3,...,则X在时刻t和y在时刻t+n之间的相关即为n阶互相关,其定义如下:

                              

            其中函数∫为计算相关系数的函数,可通过上式计算滞后n阶互相关系数的值。

            使用 airmiles和 Lakehuron数据集,来说明互相关分析的方法。其中 airmiles数据集记录了从1937年到1960年美国商业航空公司每年的乘客里程数据,而 Lakehuron数据集记录了从1875年到1972年休伦湖每年的湖平面的测量数据。将它们的时间限制在1937年到1960年,并绘制各自的时间序列曲线如图4-1所示。 

                                     

                          

                                                                                          图 4-1 

    在R语言中,可直接使用ccf函数分析时间序列的互相关性,该函数定义及参数说明如表3-1-3所示:

            

            

    plot(airmiles)
    plot(LakeHuron,xlim=c(1937,1960))
    
    ccf(airmiles, ts(LakeHuron, start = 1937, end = 1960), type="correlation")

    结果如下: 

                         

            效果如图所示,当没有延迟(即Lag=0)时,互相关系数取得最小,将近-0.8,两个时间序列大致上呈负相关关系。而当延迟10阶时,两个时间序列互相关系数最大,且为弱正相关。可通过互相关性的分析,构建用于预测的合适指标。 

    5. 典型相关分析

            典型相关是指两组变量间的相关关系,当然,它不是指对两组变量进行两两组合的简单相关,而是反映两组变量作为两个整体之间的相关性。所以典型相关的问题就是如何构建综合指标使其相关性最大。常见的做法是:

            从两组变量中提取一对线性组合,使其相关性最大;然后,从剩余的且与前面不相关的线性组合中再提取一对,使其具有最大的相关性,这样依次进行下去。所提取的线性组合就是典型变量,两组典型变量对应的相关系数就是典型相关系数。

            R语言程序包 stats中的 cancor函数,可以直接求解两个数值矩阵的典型相关系数,它的定义及参数说明如表3-1-4所示。

                  

             现使用iris数据集,分别将第1、2列和3、4列转换成数值矩阵,使用 cancor函数求解典型相关系数,代码如下。

    x<-as.matrix(iris[,1:2])
    y<-as.matrix(iris[,3:4])
    cancor(x,y)

      结果如下: 

                                     

           由结果可知,典型相关系数有两个,分别是0.9409690和0.1239369,说明第一组典型变量的相关性很强,第二组典型变量的相关性很弱。这种情况,通常第一个典型相关系数用于分析。

           xcof和ycof分别表示矩阵x和y中典型变量与原变量对应的系数矩阵,可通过标准化的数据矩阵进行乘法运算得到典型变量。

     

    展开全文
  • 在研究如何对时间序列进行线性分段的时候,浏览了60篇左右论文和教材的片段,对其中的6篇仔细阅读并编写程序和单元测试实现相应的算法。同时为了直观的看到分段效果,还制作简易的曲线图呈现原始序列和分段序列。...

     

           在研究如何对时间序列进行线性分段的时候,浏览了60篇左右论文和教材的片段,对其中的6篇仔细阅读并编写程序和单元测试实现相应的算法。同时为了直观的看到分段效果,还制作简易的曲线图呈现原始序列和分段序列。这种超负荷的工作,是在一周之内完成的,目的只有一个:选择算法。

             作为程序员,实际上并不能算是研究人员,多数情况下,他只需要不同的苹果中选择一个苹果而已,没有必要去种苹果树。

          但凡需要“选择”的时候,工作步骤如下:1、确定你想要达到的目的,这个最为重要,你的目的贯穿整个工作,千万不要在相亲的时候,突然对对方的妹妹格外关注;2、区分关注的层次比如,简要的阅读能够排除很多不需深究的东西,上面说到的60篇论文中的54篇要么是作者本身显得不妥、要么是某种方式的抄袭、要么其提供的分段图形本身就不符合要求,简单的五分钟你就能够排除,无需浪费时间。3、你感兴趣的算法各有优势和缺陷的时候,有无可能对某种主要的算法进行调整,或者组合应用其他算法的某些概念?4、实在找不到合适的算法,或者组合相应算法也无力达成的时候,能否基于你的需要而自行设计新的算法?当然,到这个层面,你也变成了那群做研究的书呆子之中的一员,不过一定要确定一点,至少你的目的明确,这和他们混稿费、混基金、呆在实验室空想是不同的,身为程序员你其实很有优势的。

            下面对算法的描述,并没有采用那些很精确的命名,而只是从算法的特征来分类。事实上大约有十来种主流的算法和近百种各类扩展、调整、优化的算法,每个都号称自己效果如何好、效率如何高、怎样支持在线划分等,但我们没有必要陷入他们的战争。选择到最后确定几种分段算法,我个人用的时间是一周,过于沉湎细节的话,恐怕一个月都无法做决断。例图中使用深圳A股深发展在2009年和2010年的实际收盘价走势,黑线为原始数据,红线为拟合线段,红点为分段点。

     

    一、对时间序列分段,是什么意思?

         时间序列,在二维平面上实际上是一条曲线,所谓分段,就是用一系列首尾相接的线段,近似的表达一条曲线。

        下图是将误差设为5级的分段图,这个粗略一些,用于长期趋势的分析。实际曲线有482个点,压缩后用22个点描绘实际曲线的趋势变化:

    误差设为5级的分段图

     

    下面是误差设为3级的分段图,和原有曲线走势更加吻合一些,用于短期波动的分析:

    误差设为3级的分段图

     

    二、为什么要对时间序列分段?

          第一个原因,很多时候,我们关注时间序列的主要趋势变化,而不太关注具体的数值和少量的异常点,对序列分段,我们可以抓住重点。比如某个公司更关心其在各城市销售数量增长没有,今天1百万,明天90万,后天110万,大后天130万,增长了没有?所以分段实际上是表达趋势,你可以将时间序列用自然语言表达,即在10月1日到11月5日之间,销售量大致保持每天增长0.5%的趋势,11月5日到11月9日之间,销售量以每日平均下跌0.02%的趋势。

         第二个原因,提高时间序列的运算性能。时间序列往往很长,数据量比较大,而序列的查询操作、聚类操作,往往涉及的计算量惊人,举例来说,2000个接近4000条记录的时间序列,遍历搜索一遍,在性能不错的笔记本上往往需要10个小时左右的时间。而缩小时间序列的长度,将大幅提高计算的效率,这个提高往往是50倍甚至100倍这样的概念。这往往让不可能的事情变得可能,比如你的一项分析需要运算一周才能出数据,届时黄花菜肯定是分外的凉,估计没有人会需要这样的软件系统。

         第三个原因,在不同精度层面上搜索。我们经常需要对同一时间序列,做长期趋势、中期趋势和短期趋势的分析。比如我们可能需要搜索60天的数据,也可能需要搜索5天的数据。通过控制划分时间序列的压缩度,我们可以在不同层次上运算,一些时候考虑长远趋势,图像会平滑一些。一些时候考虑短期趋势,图像会灵敏一些。固然,我们能够通过使用移动平均、或者将时间周期提升一个层次,来对原始序列做预处理,但在不同的压缩度的情况下划分时间序列,则是比较通用的办法。明明可以这么简单,又何须额外考虑许多数据预处理之类的工作?    

     

    三、在选择分段算法的时候,我们的核心要求是什么?

    1、可以调整搜索的精度,以获取更精细的分段和更粗略的分段。

    2、用于调整精度的阀值,要易于理解,且与不同时间序列的具体数值无关。

    3、同样的压缩比,我们希望分段的效果和原来的曲线更加吻合。

    4、次要目标:希望性能较好。这个很多时候可以放弃。

    5、次要目标:希望支持在线划分,也就是增量划分,这样我们可以保存分段的结果,新增数据后只要做少量的调整就行了。支持增量划分,则性能会达到顶点,不过要注意,前提是你要支持多少个不同精度的划分?如果算法参数过于复杂,则增量搜索要保存很多不同的分段结果,这是不合适的。所以,即使支持增量划分的算法,也要考虑你实际应用有没有可能使用的可能。

    所以对算法优劣的评估标准,应是:在相同拟合误差的情形下,压缩比越高越好;算法的阀值,不应与具体序列的数值有关。

     

    四、先说说我们结论:

    1、分段算法,可粗略的分为全局算法和局部算法。全局算法在整个时间序列的基础上,寻找最优的分段集合。局部算法根据局部的特征,从左到右寻找符合要求的分段。

    2、如果要精确表达,同时获得较好的压缩比,采用BU算法,当然这种算法的计算量较大,且不支持在线划分。如果要支持增量划分,采用SW或者层次聚类的算法。

    3、多数情况下,极值点和特征点算法,由于是局部算法,其精度和压缩比不能很好的平衡,不予考虑。

    4、分段和降噪应区分开来考虑,即先分段,然后根据需要对分段的结果降噪。这有助于简化思维,因为在分段的时候,精度和压缩比之间,本来存在矛盾,分段算法的关键在于两者间的平衡。而降噪,则是和精度、压缩比有关的,降噪会导致精度下降同时提高压缩比。

    5、阀值的确定,应该容易理解,同时也要容易选择,目前尚没有发现能够自行根据曲线的统计特征确定阀值的算法。所以我们不妨对阀值做一定的变换,比如我们将拟合误差变换成一个百分比,即误差越大则分段越少,取值在0到100之间。基于距离的,我们可以将其变换为与原始点的值的比例。算法应尽可能避免和序列的实际值有关,要保证相同的阀值针对不同序列具备相同的分段效果。

    6、永远不要考虑多维时间序列的问题,因为你能处理一维,就意味着你可以处理多维,此时,你仅需要考虑不同维度的先后秩序问题和权重问题,这个并不复杂。许多从事多维数据序列分析研究的专家,其实并不适合于做研究,因为既然涉及到这个领域并投入精力,本身说明其抽象思维方面素质太差。

     

    四、全局算法

    1、自顶向下TD算法:

        时间序列的开始点和结束点,是首先选中的分段点。然后,遍历两点之间的所有点,找出和这两点连成的直线距离最大的点,如果这个点到直线的距离“大于”预先给定的阀值,我们将其称为R,则将它作为第三个分段点。这样我们就有了两个线段,做了最初步的划分。

        之后,这个新增点到左边相邻点和右边相邻点构成的两条线段,继续寻找距离最大的点,然后,找到的两个点,谁与相应的线段距离最大,且这个距离“大于”阀值R,则该点作为第四个分段点….如此循环,直到再也找不到距离大于R的点,分段完成。

        这个阀值,也就是点到线段的距离,可以使用正交距离(原始点和分段线段在该点的值的差的绝对值)、垂直距离(原始点到分段线段的直线的长度)和欧式距离,当然也可以设置其他的特性作为阀值,比如拟合误差、又比如弧度、角度、余弦等,由此可以引申很多种不同的算法,这也是多数教授们所乐于从事的研究,简单。我们一般选择垂直距离就行了。

        这个阀值不太好理解,且与不同的时间序列具体取值有关,直接应用完全没有通用性。我本人在项目中,将其做一定的变换处理的。

     

    2、自底向上BU算法:

        这是TD算法的逆过程,首先将时间序列,划分为相邻点的短序列,当然此时的拟合误差为0,因为第一点和第二点的连线,原始点都落在线段上。将相邻两个线段连接起来,此时每条线段包含三个原始点,计算中间那个点的拟合误差。这样,所有这些三个点的线段中的中间点的拟合误差计算出来后,找出误差最小且误差小于阀值R的分段,作为第一条包含三个点的线段。

         在上面的基础上,第一条分段同样的和相邻线段连接,然后计算每一条分段的拟合误差,再找出误差最小且小于阀值R的分段,作为第二个分段。

         依次方式循环,直到所有分段的拟合误差都小于阀值R,分段结束。

         当然,你同样可以使用正交距离、垂直距离等其他属性,由此算法又演变成多种不同的算法。

     

    五、局部算法

    1、固定窗口PAA算法:

        这个最简易,阀值为R,表示你要将序列分成多少段。然后使用线段的长度除以这个值,得到窗口的长度L。从序列的第一点开始,到第L-1点作为第一段,用其平均值表示。这种算法实际上类似股票日线中的5日均价、10日均量。非常粗略,拟合不够精确。但这种算法,好象是时间序列分段研究最早的成果。

    2、滑动窗口SW算法:

       给一个窗口长度的最小值N和最大值M,然后从序列第一点开始,与第N点连线,计算各点拟合误差,如果误差总值小于给定的阀值R,拟合成功,增加窗口的长度,连线再计算拟合误差,如果拟合误差小于R,继续增加窗口长度一直到窗口长度为M。如果拟合误差大于R,则第一段结束,该点作为新的窗口的起始点,继续同样的过程,直到序列划分完毕。

          这种算法属于局部算法,有方向性,从算法的原理来看,就容易漏报和将本应拟合的线段划分为两条线段。代码实现的效果来看,划分的不太理想。

    3、极值点:

        所谓极值点,是指曲线中的转折点,寻找极值点的算法很简单,从序列第二点开始遍历,没点如果同时比前后两点大,或者同时比前后两点小,那么这个点就是时间序列的一个极值点。找出所有极值点的集合,就得到了所有分段点。这种算法,是基于人类视角的关注点而来,但很显然存在两个问题:1、是局部算法,只考虑与相邻两点的问题,所以划分的时候容易过于精确而忽略趋势。2、只考虑了趋势反向问题,即从上升到下降和从下降到上升,没有考虑趋势加速问题,比如先是15%角度上升,然后是40%角度上升,之间的结合点显然漏掉了。这也意味着此种算法从原理上就不能精确表达趋势。优点也是有的,即计算简单,性能很好,支持在线划分,没有阀值。

    4、特征点:

        在极值点的基础上,加上判断,比如两个极值点之间的距离必须大于某个阀值,或者该极值点除了比前后点同时大或同时小之外,还需要大于某个比例阀值R。这样的改进意义不大,只是进一步减少极值点而已。

    5、趋势点:

        在极值点的基础上,考虑了趋势加速问题,即计算某点和相邻点的角度,或者弧度、或者余弦等,大于某个特定阀值R的称为趋势点,按照这种算法基本上不会遗漏掉极值点,趋势加速或减速的关键点也不会遗漏掉。不过这同样是一种局部算法,极值点方法其他毛病,这种算法并没有避免。

    6、层次聚类:实际上仅仅运用了KMean算法的中心和半径的概念,其他没多少区别。这种算法的阀值比较多,而且看起来分类的效果不好,没有看点。不过这种算法,其分类的效果好于SW算法,也同样的支持在线划分。

     

    六、时间序列自身特征问题:

        一部分类型的时间序列,遵循一定的规律,价格逐日波动的过程相对平稳,另一些很可能出现前后两天变动数倍的情形。

        这样,首先我们面临的问题,是能否针对不同业务含义的数据调整精度,以正常的分段;第二,是能否对数据进行预处理,以支持分段。

        这两种方式,事实上都与业务系统相关。

        反过来说,我们需要研究,能否通过自动分析数据本身的特征,由算法来决定是否需要预处理、基于何种模式设置何种精度合适。这里涉及到的序列特征,包括序列的极值(最大与最小的差距)、序列的标准差(反映离散程度)、序列的最大取值范围(出现最多的区域,出现百分比多少),当然,这仅仅是一些统计学方面的概念。

        事实上,对于波幅巨大的序列,如果分段,则很可能出现全部原始点都是分段点的情形,这种情形下的趋势价值有限。使用移动平均的方式来平滑是一个办法,但这里的相似性度量仍然会存在问题。

     

    七、进一步关注的内容:

        时间序列方面的具体应用,无非落实在识别、趋势分析和异常检测,在趋势分析方面,基于拟合的思路远不如使用时间序列基于搜索、聚类和统计的思路。基于拟合的思路,做决策支持分析的时候,是无法说服客户的。比如神经网络之类的方式,本质上就是用公式表示现在的走势,同时认为未来的走势也符合这个公式的规律,这个显然是期望使用数学来统治世界,很难令人信服。其他诸如ARMA之类的方法,大体上也是同样的路子。 只有基于搜索、聚类和统计的方法,才能够有一定的说服力,因为这里找出的是一些历史上已经发现的模式,期间蕴藏着我们不容易表达的规律。

          由于时间序列,或者说有向序列,通常的需求无非在趋势和异常两点。而表达趋势方面,需要一定的模糊性,那么模糊数学、机器学习方面的一些概念,能否用于这个领域,目前尚没有特别有效的方法。

          现有的时间序列工具,当图像被当作二维平面上的点的时候,按照先左后右获得第一条序列,从上到下获得所有序列这样的顺序,是能够进行分析的,比如找出相片中有多少人、找出相片中的人是不是你本人、找出相片中出现了一辆车没有、指纹是否符合等、雷达中获取的反射信息构成的模糊形状是否可能是一架飞机?当然,同样的扩展到三维空间,则对3d图形实际上也能够处理,无非是确定序列数量和方向、确定各序列权重的问题。

    转载于:https://www.cnblogs.com/by1990/archive/2011/01/15/1936296.html

    展开全文
  • 时间序列分析

    千次阅读 2018-10-29 21:15:51
    时间序列 时间序列(简称为时序)是指...时间序列分析是一种动态数据处理的统计方法。该方法基于随机过程理论和数理统计学方法,研究随机数据序列所遵从的统计变化规律,以用于解决实际问题。通常影响时间序列变化的...
  • 时间序列

    千次阅读 2018-09-02 00:33:43
    时间序列进行观察、研究、找寻他发展变化的规律,预测他将来的走势就是时间序列分析,时间序列分析方法只适用于近期与短期的预测。 相关特征统计量: 均值函数序列:反映的是时间序列每时每刻的平均水平 方差...
  • 时间序列预测方法最全总结!

    千次阅读 多人点赞 2021-03-12 00:15:38
    时间序列预测就是利用过去一段时间的数据来预测未来一段时间内的信息,包括连续型预测(数值预测,范围估计)与离散型预测(事件预测)等,具有非常高的商业价值。需要明确一点的是,与回归分析预测模型...
  • 线性插补的方法包括一元线性回归,多元线性回归,岭回归等,有关方法可以参考之前的文章:《TS技术课堂 | 时间序列回归》 季节性+线性插值 经济数据或者季节波动数据,常常不符合简单的线性变化强假设。对这样的数据...
  • 时间序列预测方法总结

    千次阅读 2020-10-18 15:14:52
    这本来是我回答的一个问题:有什么好的模型可以做高精度的时间序列预测呢? - BINGO Hong的回答 - 知乎 https://www.zhihu.com/question/21229371/answer/533770345 但觉得在那个答案下一直更新好麻烦,干脆就移到...
  • 时间序列预测算法总结

    万次阅读 多人点赞 2018-10-18 10:30:48
    时间序列算法 time series data mining 主要包括decompose(分析数据的各个成分,例如趋势,周期性),prediction(预测未来的值),classification(对有序数据序列的feature提取与分类),clustering(相似数列...
  • 时间序列分析方法——ARIMA模型案例

    千次阅读 2020-06-19 18:35:24
    目录一、方法简介数据示例二、ARIMA模型python建模过程[^2]1 添加基础库2 读取数据3 绘制时间序列图4 自相关5 平稳性检验6 时间序列的差分d7 合适的p,q8 模型检验Ljung-Box检验9 模型预测 时间序列分析方法1主要有...
  • 时间序列分析之相关性

    万次阅读 2019-02-20 17:27:27
    两种时间序列的相关性 方差 (Variance) 设随机变量X的均值 E(X) = m,则描述 X 的取值和它的均值 m 之间的偏差程度大小的数字特征就是方差。 但是不能直接用 E(X - m) 来表示方差,因为 E...
  • 时间序列相关函数

    千次阅读 2015-08-25 17:16:13
    文档1:《R与金融时间序列分析常见问题集》 【包】 library(zoo) #时间格式预处理 library(xts) #同上 library(timeSeires) #同上 library(urca) #进行单位根检验 library(tseries) #arma模型 lib
  • 时间序列matlab的实现

    万次阅读 多人点赞 2019-08-23 11:14:45
    时间序列预测法其实是一种回归预测方法,属于定量预测,其基本原理是:一方面承认事物发展的延续性,运用过去的时间序列数据进行统计分析,推测出事物的发展趋势;另一方面充分考虑到由于偶然因素影响而产生的随机性...
  • 时间序列分析之ADF检验

    万次阅读 多人点赞 2019-02-06 18:59:09
    在使用很多时间序列模型的时候,如 ARMA、ARIMA,都会要求时间序列是平稳的,所以一般在研究一段时间序列的时候,第一步都需要进行平稳性检验,除了用肉眼检测的方法,另外比较常用的严格的统计检验方法就是ADF检验...
  • 随着时间推移,制造设备比如贴片机的位置由于各种原因会产生小的偏差。这些偏差可能是阶跃,也有可能是渐变的形式。由于偏差值很小,产线的自动光学检测设备并不会报警;然而小的偏差如果不经处理,经过一定时间累积...
  • 时间序列, 智能运维时间序列的自回归模型—从线性代数的角度来看MARCH 23, 2018 LEAVE A COMMENTFibonacci 序列在开始讲解时间序列的自回归模型(AutoRegression Model)之前,我们需要先介绍一下线性代数的基础...
  • 时间序列预测,非季节性ARIMA及季节性SARIMA

    万次阅读 多人点赞 2019-03-24 21:55:00
    我们将首先介绍和讨论自相关,平稳性和季节性的概念,并继续应用最常用的时间序列预测方法之一,称为ARIMA。 介绍 时间序列提供了预测未来价值的机会。基于以前的价值观,可以使用时间序列来预测经济,天气和...
  • 时间序列的平稳性、随机性检验 在拿到时间序列数据后,首先要对数据的随机性和平稳性进行检测, 这两个检测是时间序列预测的重要部分。根据不同检测结果需要采取不同的分析方法。 为什么时间序列要求平稳性呢?平稳...
  • 数据挖掘之时间序列分析

    万次阅读 多人点赞 2018-08-12 23:45:16
    按时间顺序排列的一组随机变量X1,X2,…,Xt表示一个随机事件的时间序列时间序列分析的目的是给定一个已被观测了的时间序列,预测该序列的未来值。 表1 常用的时间序列模型 模型名称 描述 平滑法 常...
  • 最近一直在接触时间序列,所以打算写一些有关时间序列的文章,预测部分会从规则开始、到传统模型、到机器学习、再到深度学习,此外也会介绍一些时间序列的基本概念,包括自相关、平稳性、滞后性、季节...
  • 时间序列预测

    千次阅读 多人点赞 2021-01-25 15:20:33
    1、时间序列是跟时间有关的,而线性回归模型中观察结果是独立的; 2、随着上升或者下降的趋势,更多的时间序列出现季节性趋势的形式。 常用按时间序列排列的一组随机变量X_1,X_2,…,X_t来表示一个随机事件序列,记为...
  • 13、时间序列

    千次阅读 2016-02-12 10:56:50
    1、时间序列的分类 时间序列分为:非季节性数据和季节性数据 一个非季节性时间序列包含一个趋势部分和一个不规则部分。 一个季节性时间序列包含一个趋势部分,一个季节性部分和一个不规则部分。 在实践操作中,时间...
  • R时间序列分析

    千次阅读 2016-06-10 21:05:17
    R时间序列分析 为什么定阶数,如何定,如何判断 R时间序列分析工具 xts包 xts(x=NUll,order.by=index(x),…) coredata() xts数据子集 OHLC数据格式 quantmod包 TTR包 自回归模型(AR) 跟以前时刻有关和当前随机...
  • 时间序列分析之协整检验

    万次阅读 多人点赞 2019-02-07 14:18:02
    平稳性是进行时间序列分析的一个很重要的前提,很多模型都是基于平稳下进行的,而现实中,很多时间序列都是非平稳的,所以协整是从分析时间序列的非平稳性入手的。 协整的内容是: 设序列是 d 阶单整的,记为,...
  • 如何判断时间序列属于加法模型还是乘法模型 如果时间序列图的趋势随着时间的推移,序列的季节波动变得越来越大,则建议使用乘法模型;如果序列的季节波动能够基本维持恒定,则建议使用加法模型。 时间序列的预测步骤...
  • 时间序列初级理论篇

    千次阅读 2017-10-23 09:11:52
    前言 数学特征 1 一般随机变量的数学特征 11 期望 12 方差 13 协方差 ...111 弱宽平稳时间序列的数学特征 2 差分法 3 其他变换 时间序列影响因素影响因素的叠加 平稳时间序列模型介绍 1 自回
  • R语言实现时间序列分析

    千次阅读 2019-11-11 20:11:45
    一、时间序列分析导论图: 处理方式: 1.对于原始数据进行季节性处理和差分,以形成平稳序列;期间如果遇到了随机序列,则停止时间序列建模 2.对于给定的序列进行自相关函数和偏自相关函数分析(在不同的滞后...

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 29,798
精华内容 11,919
关键字:

判断时间序列线性相关的方法