time-series - 在 Stata / SAS 中创建具有缺失值的运行平均值

Question

多年来，我对环境和气象变量（温度和湿度）进行了每小时测量的时间序列。从这些每小时值中，我想计算 24 小时运行平均值来创建曝光参数。为此，要求至少有 17 个每小时测量值可用，且连续缺失值不超过 6 小时。如果 24 小时内连续缺失超过 6 个小时值，则该特定日期的数据将设置为缺失。如何在 Stata 或 SAS 中实现这一点？

提前致谢

Answer 1

看起来您可以使用以下组合为“有效”观察创建一个虚拟变量

• by varname : generate ....,

• egen， 和

• 滞后变量 ( L.varname, L2.varname... L24.varname...)

然后，使用您的数据子集创建您的平均值（例如， yourcommand ... if dummy==1 ...）

Answer 2

一个Stata解决方案：

1. 用于tssmooth ma myvar_ma = myvar, w(24)创建移动平均线。缺失的将被忽略。

2. 创建指标gen ismiss = missing(myvar)

3. 用于tssmooth ma ismiss_ma = ismiss, w(24)创建指标的移动平均线。

4. replace myvar_ma = . if ismiss_ma > (7/24)

（至少 17/24 必须存在，因此 7 或更少的缺失是可以接受的，但 8 或更多是不可接受的。

编辑。tsegenSSC 现在提供了一种解决此类问题的简单方法。您可以直接在命令语法中指定窗口中可接受的最小非缺失值数量。

Answer 3

好的，这是我的尝试。首先创建一些示例数据以供使用：

**
** CREATE ~3 YEARS DAYS OF HOURLY TEMPERATURE DATA
** THIS IS UGLY - IM SURE THERES A BETTER WAY TO DO IT BUT WHATEVER
*;
data tmp;
  pi = constant('pi');
  do year=1 to 3;
    linear_trend = 0.1 * year;
    day = 0;
    do yearly_progress=0 to (pi*2) by (pi/182.5);
      day = day + 1;
      yearly_seasonality = (1 + sin(yearly_progress)) / 2;
      hour = 0;
      day_mod = (ranuni(0)*10);
      do hourly_progress=0 to (pi*2) by (pi/12);
        hourly_seasonality = (1 + sin(hourly_progress)) / 2;
        if hour ne 24 then do;
          temperature = 60*(1+linear_trend) + (20 * yearly_seasonality) + (30 * hourly_seasonality) - day_mod;
          output;
        end;
        hour = hour + 1;
      end;
    end;
  end;
run;


**
** ADD SOME MISSING VALS
** ~ 10% MISSING
** ~ 10 IN A ROW MISSING EVERY 700 OR SO HOURS
*;
data sample_data;
  set tmp;
  if (ranuni(0) < 0.1) or (mod(_n_,710) > 700) then do;
    temperature = .;
  end;
run;

其次，如果满足要求，则计算温度的移动平均值：

**
** I DONT NORMALLY LIKE USING THE LAG FUNCTION BUT IN THIS CASE ITS IDEAL
*;
data results;
  set sample_data;

  **
  ** POPULATE AN ARRAY WITH THE 24 CURRENT VALUES
  ** BECAUSE WE ARE USING LAG FUNCTION MAKE SURE IT IS NOT WITHIN ANY 
  ** CONDITIONAL IF STATEMENTS
  *;
  array arr [0:23] temperature0-temperature23;
  temperature0  =  lag0(temperature);
  temperature1  =  lag1(temperature);
  temperature2  =  lag2(temperature);
  temperature3  =  lag3(temperature);
  temperature4  =  lag4(temperature);
  temperature5  =  lag5(temperature);
  temperature6  =  lag6(temperature);
  temperature7  =  lag7(temperature);
  temperature8  =  lag8(temperature);
  temperature9  =  lag9(temperature);
  temperature10 = lag10(temperature);
  temperature11 = lag11(temperature);
  temperature12 = lag12(temperature);
  temperature13 = lag13(temperature);
  temperature14 = lag14(temperature);
  temperature15 = lag15(temperature);
  temperature16 = lag16(temperature);
  temperature17 = lag17(temperature);
  temperature18 = lag18(temperature);
  temperature19 = lag19(temperature);
  temperature20 = lag20(temperature);
  temperature21 = lag21(temperature);
  temperature22 = lag22(temperature);
  temperature23 = lag23(temperature);

  **
  ** ITERATE OVER THE ARRAY VARIABLES TO MAKE SURE WE MEET THE REQUIREMENTS
  *;
  available_observations  = 0;
  missing_observations    = 0;
  max_consecutive_missing = 0;
  tmp_consecutive_missing = 0;
  do i=0 to 23;
    if arr[i] eq . then do;
      missing_observations    = missing_observations + 1;
      tmp_consecutive_missing = tmp_consecutive_missing + 1;
      max_consecutive_missing = max(max_consecutive_missing, tmp_consecutive_missing);
    end;
    else do;
      available_observations  = available_observations + 1;        
      tmp_consecutive_missing = 0;
    end;
  end;

  if tmp_consecutive_missing <= 6 and available_observations >= 17 then do;
    moving_avg = mean(of temperature0-temperature23);
  end;
run;

Answer 4

对于一般移动平均计算，使用 PROC EXPAND 是最简单的方法（您需要获得 ETS 许可才能使用此程序）。例如，下面的代码将计算 24 个周期的移动平均值并将前 16 个观测值设置为缺失。但是，为了符合您的其他标准，您之后仍需要按照 Rob 的代码行运行数据步骤，因此您不妨在该步骤中执行所有计算。

proc expand data=sample_data out=mov_avg_results;
convert temperature=mean_temp / method=none transformout=(movave 24 trimleft 17);
run;