估計階段和季節週期

Question

，我有以下數據的幅度：

CET <- url("http://www.metoffice.gov.uk/hadobs/hadcet/cetml1659on.dat") 
cet <- read.table(CET, sep = "", skip = 6, header = TRUE, 
        fill = TRUE, na.string = c(-99.99, -99.9)) 
names(cet) <- c(month.abb, "Annual") 
cet <- cet[-nrow(cet), ] 
rn <- as.numeric(rownames(cet)) 
Years <- rn[1]:rn[length(rn)] 
annCET <- data.frame(Temperature = cet[, ncol(cet)],Year = Years) 
cet <- cet[, -ncol(cet)] 
cet <- stack(cet)[,2:1] 
names(cet) <- c("Month","Temperature") 
cet <- transform(cet, Year = (Year <- rep(Years, times = 12)), 
       nMonth = rep(1:12, each = length(Years)), 
       Date = as.Date(paste(Year, Month, "15", sep = "-"),format = "%Y-%b-%d")) 
cet <- cet[with(cet, order(Date)), ] 

idx <- cet$Year > 1900 
cet <- cet[idx,] 
cet <- cet[,c('Date','Temperature')] 

plot(cet, type = 'l') 

這表明從1900年每月的溫度循環至2014年在英國，英國。

我想通過this論文中概述的方法來評估溫度季節週期的相位和振幅。具體而言，它們描述了給定的12個月的值（如我們在這裏），我們可以估算每年成分爲：

，其中X（t）表示表面溫度的12個月的值，X（T + T0 ），t = 0.5，...，11.5，是每月的12個月的溫度值，其中2的因子是考慮正頻率和負頻率。

然後季節性週期的振幅和相位可以計算爲

和

他們指定，這些數據每年，他們計算出年度（每年一個週期）使用傅立葉變換的正弦分量，如以上所示的等式。

我有點卡在如何生成他們在這裏演示的時間序列。任何人都可以提供一些關於如何重現這些方法的指導。請注意，我也在matlab中工作 - 以防萬一任何人對如何在該環境中實現這一點有一些建議。

這是數據的一個子集。

Date Temperature 
1980-01-15 2.3 
1980-02-15 5.7 
1980-03-15 4.7 
1980-04-15 8.8 
1980-05-15 11.2 
1980-06-15 13.8 
1980-07-15 14.7 
1980-08-15 15.9 
1980-09-15 14.7 
1980-10-15 9 
1980-11-15 6.6 
1980-12-15 5.6 
1981-01-15 4.9 
1981-02-15 3 
1981-03-15 7.9 
1981-04-15 7.8 
1981-05-15 11.2 
1981-06-15 13.2 
1981-07-15 15.5 
1981-08-15 16.2 
1981-09-15 14.5 
1981-10-15 8.6 
1981-11-15 7.8 
1981-12-15 0.3 
1982-01-15 2.6 
1982-02-15 4.8 
1982-03-15 6.1 
1982-04-15 8.6 
1982-05-15 11.6 
1982-06-15 15.5 
1982-07-15 16.5 
1982-08-15 15.7 
1982-09-15 14.2 
1982-10-15 10.1 
1982-11-15 8 
1982-12-15 4.4 
1983-01-15 6.7 
1983-02-15 1.7 
1983-03-15 6.4 
1983-04-15 6.8 
1983-05-15 10.3 
1983-06-15 14.4 
1983-07-15 19.5 
1983-08-15 17.3 
1983-09-15 13.7 
1983-10-15 10.5 
1983-11-15 7.5 
1983-12-15 5.6 
1984-01-15 3.8 
1984-02-15 3.3 
1984-03-15 4.7 
1984-04-15 8.1 
1984-05-15 9.9 
1984-06-15 14.5 
1984-07-15 16.9 
1984-08-15 17.6 
1984-09-15 13.7 
1984-10-15 11.1 
1984-11-15 8 
1984-12-15 5.2 
1985-01-15 0.8 
1985-02-15 2.1 
1985-03-15 4.7 
1985-04-15 8.3 
1985-05-15 10.9 
1985-06-15 12.7 
1985-07-15 16.2 
1985-08-15 14.6 
1985-09-15 14.6 
1985-10-15 11 
1985-11-15 4.1 
1985-12-15 6.3 
1986-01-15 3.5 
1986-02-15 -1.1 
1986-03-15 4.9 
1986-04-15 5.8 
1986-05-15 11.1 
1986-06-15 14.8 
1986-07-15 15.9 
1986-08-15 13.7 
1986-09-15 11.3 
1986-10-15 11 
1986-11-15 7.8 
1986-12-15 6.2 
1987-01-15 0.8 
1987-02-15 3.6 
1987-03-15 4.1 
1987-04-15 10.3 
1987-05-15 10.1 
1987-06-15 12.8 
1987-07-15 15.9 
1987-08-15 15.6 
1987-09-15 13.6 
1987-10-15 9.7 
1987-11-15 6.5 
1987-12-15 5.6 
1988-01-15 5.3 
1988-02-15 4.9 
1988-03-15 6.4 
1988-04-15 8.2 
1988-05-15 11.9 
1988-06-15 14.4 
1988-07-15 14.7 
1988-08-15 15.2 
1988-09-15 13.2 
1988-10-15 10.4 
1988-11-15 5.2 
1988-12-15 7.5 
1989-01-15 6.1 
1989-02-15 5.9 
1989-03-15 7.5 
1989-04-15 6.6 
1989-05-15 13 
1989-06-15 14.6 
1989-07-15 18.2 
1989-08-15 16.6 
1989-09-15 14.7 
1989-10-15 11.7 
1989-11-15 6.2 
1989-12-15 4.9 
1990-01-15 6.5 
1990-02-15 7.3 
1990-03-15 8.3 
1990-04-15 8 
1990-05-15 12.6 
1990-06-15 13.6 
1990-07-15 16.9 
1990-08-15 18 
1990-09-15 13.2 
1990-10-15 11.9 
1990-11-15 6.9 
1990-12-15 4.3 
1991-01-15 3.3 
1991-02-15 1.5 
1991-03-15 7.9 
1991-04-15 7.9 
1991-05-15 10.8 
1991-06-15 12.1 
1991-07-15 17.3 
1991-08-15 17.1 
1991-09-15 14.7 
1991-10-15 10.2 
1991-11-15 6.8 
1991-12-15 4.7 
1992-01-15 3.7 
1992-02-15 5.4 
1992-03-15 7.5 
1992-04-15 8.7 
1992-05-15 13.6 
1992-06-15 15.7 
1992-07-15 16.2 
1992-08-15 15.3 
1992-09-15 13.4 
1992-10-15 7.8 
1992-11-15 7.4 
1992-12-15 3.6 
1993-01-15 5.9 
1993-02-15 4.6 
1993-03-15 6.7 
1993-04-15 9.5 
1993-05-15 11.4 
1993-06-15 15 
1993-07-15 15.2 
1993-08-15 14.6 
1993-09-15 12.4 
1993-10-15 8.5 
1993-11-15 4.6 
1993-12-15 5.5 
1994-01-15 5.3 
1994-02-15 3.2 
1994-03-15 7.7 
1994-04-15 8.1 
1994-05-15 10.7 
1994-06-15 14.5 
1994-07-15 18 
1994-08-15 16 
1994-09-15 12.7 
1994-10-15 10.2 
1994-11-15 10.1 
1994-12-15 6.4 
1995-01-15 4.8 
1995-02-15 6.5 
1995-03-15 5.6 
1995-04-15 9.1 
1995-05-15 11.6 
1995-06-15 14.3 
1995-07-15 18.6 
1995-08-15 19.2 
1995-09-15 13.7 
1995-10-15 12.9 
1995-11-15 7.7 
1995-12-15 2.3 
1996-01-15 4.3 
1996-02-15 2.5 
1996-03-15 4.5 
1996-04-15 8.5 
1996-05-15 9.1 
1996-06-15 14.4 
1996-07-15 16.5 
1996-08-15 16.5 
1996-09-15 13.6 
1996-10-15 11.7 
1996-11-15 5.9 
1996-12-15 2.9 
1997-01-15 2.5 
1997-02-15 6.7 
1997-03-15 8.4 
1997-04-15 9 
1997-05-15 11.5 
1997-06-15 14.1 
1997-07-15 16.7 
1997-08-15 18.9 
1997-09-15 14.2 
1997-10-15 10.2 
1997-11-15 8.4 
1997-12-15 5.8 
1998-01-15 5.2 
1998-02-15 7.3 
1998-03-15 7.9 
1998-04-15 7.7 
1998-05-15 13.1 
1998-06-15 14.2 
1998-07-15 15.5 
1998-08-15 15.9 
1998-09-15 14.9 
1998-10-15 10.6 
1998-11-15 6.2 
1998-12-15 5.5 
1999-01-15 5.5 
1999-02-15 5.3 
1999-03-15 7.4 
1999-04-15 9.4 
1999-05-15 12.9 
1999-06-15 13.9 
1999-07-15 17.7 
1999-08-15 16.1 
1999-09-15 15.6 
1999-10-15 10.7 
1999-11-15 7.9 
1999-12-15 5 
2000-01-15 4.9 
2000-02-15 6.3 
2000-03-15 7.6 
2000-04-15 7.8 
2000-05-15 12.1 
2000-06-15 15.1 
2000-07-15 15.5 
2000-08-15 16.6 
2000-09-15 14.7 
2000-10-15 10.3 
2000-11-15 7 
2000-12-15 5.8 
2001-01-15 3.2 
2001-02-15 4.4 
2001-03-15 5.2 
2001-04-15 7.7 
2001-05-15 12.6 
2001-06-15 14.3 
2001-07-15 17.2 
2001-08-15 16.8 
2001-09-15 13.4 
2001-10-15 13.3 
2001-11-15 7.5 
2001-12-15 3.6 
2002-01-15 5.5 
2002-02-15 7 
2002-03-15 7.6 
2002-04-15 9.3 
2002-05-15 11.8 
2002-06-15 14.4 
2002-07-15 16 
2002-08-15 17 
2002-09-15 14.4 
2002-10-15 10.1 
2002-11-15 8.5 
2002-12-15 5.7 
2003-01-15 4.5 
2003-02-15 3.9 
2003-03-15 7.5 
2003-04-15 9.6 
2003-05-15 12.1 
2003-06-15 16.1 
2003-07-15 17.6 
2003-08-15 18.3 
2003-09-15 14.3 
2003-10-15 9.2 
2003-11-15 8.1 
2003-12-15 4.8 
2004-01-15 5.2 
2004-02-15 5.4 
2004-03-15 6.5 
2004-04-15 9.4 
2004-05-15 12.1 
2004-06-15 15.3 
2004-07-15 15.8 
2004-08-15 17.6 
2004-09-15 14.9 
2004-10-15 10.5 
2004-11-15 7.7 
2004-12-15 5.4 
2005-01-15 6 
2005-02-15 4.3 
2005-03-15 7.2 
2005-04-15 8.9 
2005-05-15 11.4 
2005-06-15 15.5 
2005-07-15 16.9 
2005-08-15 16.2 
2005-09-15 15.2 
2005-10-15 13.1 
2005-11-15 6.2 
2005-12-15 4.4 
2006-01-15 4.3 
2006-02-15 3.7 
2006-03-15 4.9 
2006-04-15 8.6 
2006-05-15 12.3 
2006-06-15 15.9 
2006-07-15 19.7 
2006-08-15 16.1 
2006-09-15 16.8 
2006-10-15 13 
2006-11-15 8.1 
2006-12-15 6.5 
2007-01-15 7 
2007-02-15 5.8 
2007-03-15 7.2 
2007-04-15 11.2 
2007-05-15 11.9 
2007-06-15 15.1 
2007-07-15 15.2 
2007-08-15 15.4 
2007-09-15 13.8 
2007-10-15 10.9 
2007-11-15 7.3 
2007-12-15 4.9 
2008-01-15 6.6 
2008-02-15 5.4 
2008-03-15 6.1 
2008-04-15 7.9 
2008-05-15 13.4 
2008-06-15 13.9 
2008-07-15 16.2 
2008-08-15 16.2 
2008-09-15 13.5 
2008-10-15 9.7 
2008-11-15 7 
2008-12-15 3.5 
2009-01-15 3 
2009-02-15 4.1 
2009-03-15 7 
2009-04-15 10 
2009-05-15 12.1 
2009-06-15 14.8 
2009-07-15 16.1 
2009-08-15 16.6 
2009-09-15 14.2 
2009-10-15 11.6 
2009-11-15 8.7 
2009-12-15 3.1 
2010-01-15 1.4 
2010-02-15 2.8 
2010-03-15 6.1 
2010-04-15 8.8 
2010-05-15 10.7 
2010-06-15 15.2 
2010-07-15 17.1 
2010-08-15 15.3 
2010-09-15 13.8 
2010-10-15 10.3 
2010-11-15 5.2 
2010-12-15 -0.7 
2011-01-15 3.7 
2011-02-15 6.4 
2011-03-15 6.7 
2011-04-15 11.8 
2011-05-15 12.2 
2011-06-15 13.8 
2011-07-15 15.2 
2011-08-15 15.4 
2011-09-15 15.1 
2011-10-15 12.6 
2011-11-15 9.6 
2011-12-15 6 
2012-01-15 5.4 
2012-02-15 3.8 
2012-03-15 8.3 
2012-04-15 7.2 
2012-05-15 11.7 
2012-06-15 13.5 
2012-07-15 15.5 
2012-08-15 16.6 
2012-09-15 13 
2012-10-15 9.7 
2012-11-15 6.8 
2012-12-15 4.8 
2013-01-15 3.5 
2013-02-15 3.2 
2013-03-15 2.7 
2013-04-15 7.5 
2013-05-15 10.4 
2013-06-15 13.6 
2013-07-15 18.3 
2013-08-15 16.9 
2013-09-15 13.7 
2013-10-15 12.5 
2013-11-15 6.2 
2013-12-15 6.3 
2014-01-15 5.7 
2014-02-15 6.2 
2014-03-15 7.6 
2014-04-15 10.2 
2014-05-15 12.2 
2014-06-15 15.1 
2014-07-15 17.7 
2014-08-15 14.9 
2014-09-15 15.1 
2014-10-15 12.5 
2014-11-15 8.6 
2014-12-15 5.2 

Answer 1

從字面上看，爲Y公式可以在MATLAB中表示爲：

t=0.5:0.5:11.5; %//make sure the step size is indeed 0.5 
Y = 1/6.*sum(exp(2*pi*i.*t/12).*X(t0-t); %// add the function for X 
phi = atan2(imag(Y)/real(Y)); %// seasonal phase 

不知道爲X功能我不能肯定這的確可以矢量化，或者你是否會有循環，可以像做：在你想要的任何t0

t=0.5:0.5:11.5; %//make sure the step size is indeed 0.5 
Ytmp(numel(t),1)=0; %// initialise output 
for ii = 1:numel(t) 
    Ytmp(ii,1) = exp(2*pi*i.*t(ii)/12).*X(t0-t(ii)); 
end 
Y = 1/6 * sum(Ytmp) 

只是插槽，環比上面的代碼，你有你的時間序列。