How to calculate trend value per year for time series data?

I have data of water level (both maximum and minimum in a scale of the meter) on a daily basis for a certain number of years. The data set has some missing values as well. I have replaced those missing values by NaN. I want to calculate the trend in changes of water level per year. I also need to compute the standard error, to see if the trends are statistically significant or not. Is there any MATLAB functions available for doing this analysis? I used the 12 months running mean to calculate the annual mean for removing the seasonal cycles (codes are as follows).

alldays=dados.textdata; 
[ntotal nada]=size(alldays)
water=dados.data; 
[nvalues nada]=size(water)
max_monthly=water(:,1);
min_monthly=water(:,2);
bissexto=[1952 1956 1960 1964 1968 1972 1976 1980 1984 1988 1992 1996 2000 2004 2008 2012 2016 2020 2024 2028 2032 2036 2040 2044 2048 2052 2056 2060 2064 2068 2072 2076 2080 2084 2088 2092 2096];

figure
plot(max_monthly,'r-')
hold on
plot(min_monthly,'b-')
ylim([-4 6])

% close all
%%%%%%%%Starting day %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
day=1;
month=4;
year=1990;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
m=1;
k=1;
p=1;

nyear=1;
[tam nada]=size(max_monthly);
d=1;
som1=max_monthly(1);
som2=min_monthly(1);
contador=2;
while (d<tam)
  d=d+1; 
  f=find (year==bissexto);
  TF=isempty(f);
  if(TF==1)
    ndays=[31 28 31 30 31 30 31 31 30 31 30 31];
  else
    ndays=[31 29 31 30 31 30 31 31 30 31 30 31];
    %         ano
  end
  if(day>=ndays(month))
    soma1=nansum(som1);
    soma2=nansum(som2);
    avg1=soma1./ndays(month);
    avg2=soma2./ndays(month);
    %            if avg1==0
    %                avg1=NaN;
    %            end
    %            if avg2==0
    %                avg2=NaN;
    %            end
    if(month<12)
      max_monthly(m)=[avg1];
      min_monthly(m)=[avg2];
      max_allmonths(p)=avg1;
      min_allmonths(p)=avg2;
      clear avg1 avg2 contador som1 som2 soma1 soma2
      contador=1;
      m=m+1;
      k=k+1;
      p=p+1;
    end
    day=1;
    som1=max_monthly(d);
    som2=min_monthly(d);
    contador=contador+1;
    month=month+1;
  else
    som1(contador)=max_monthly(d);
    som2(contador)=min_monthly(d);
    contador=contador+1;
    day=day+1;
  end
  if(month>12)
    max_year(nyear)=mean(max_allmonths);
    min_year(nyear)=mean(min_allmonths);
    k=1;
    month=1;
    year=year+1;
    nyear=nyear+1;
  else
    month=month;
    year=year;
  end
  clear ndays TF f figura3
 end
 %%Difference
 DifMaxMin=max_year-abs(min_year);
figure
plot(max_year,'r-')
hold on
plot(min_year,'b-')
figure
subplot(2,1,1)
plot(1:tam,max_monthly,'r-')
hold on
plot(min_monthly,'b-')
subplot(2,1,2)
plot(1:length(max_year),max_year,'r-')
hold on
plot(1:length(min_year),min_year,'b-')

I am very new in MATLAB. So, any kind of help will be appreciated immensely! Thanks again!

[dpb]

Import into a timetable and use retime and you'll be done in just a line or two...

[Bishwajit Roy]

Hi, I revised my question. Can you be bit specific or provide me the functions or code, please? Thanks! Roy

Answer 1

[dpb]

For starters using timetable and retime...

opts=detectImportOptions('barguna.txt');               % build basic import options object
opts=setvartype(opts,1,'datetime');                    % set the first column datetime type
opts=setvaropts(opts,1,'InputFormat','dd-MM-yyyy');    % and the input format for it
opts.VariableNames={'Date','MaxLev','MinLev'};         % define useful variable names
tt=table2timetable(readtable('barguna.txt',opts));     % and read table, convert to timetable
yrmeans=retime(tt,'yearly',@nanmean);                  % compute the annual means
yr=yrmeans.Date; yr.Format='yyyy'; yrmeans.Date=yr;    % change to only display year

>> yrmeans                                             % see what we got...
yrmeans =
29×2 timetable
  Date    MaxLev     MinLev  
  ____    ______    _________
  1990    2.0609     0.090364
  1991    1.8643    -0.094959
  1992     1.876    -0.056967
  1993    1.9142    -0.013014
  1994    1.8547     -0.12948
  1995    1.8788    -0.066192
  1996    1.8744    -0.076093
  1997    1.7816     -0.18792
  1998    1.8968     0.025753
  1999    2.0535    -0.046767
  2000    2.0557    -0.028197
  2001    2.1407      0.11874
  2002    2.0523    -0.097019
  2003    2.3064      0.18219
  2004    1.6076      0.30488
  2005    1.6687      0.53648
  2006    1.6655      0.36096
  2007     1.805      0.31266
  2008    1.8197      0.01776
  2009    1.9946     0.025644
  2010    2.0095      0.17868
  2011    2.0188     0.096154
  2012    1.8374      0.15153
  2013    2.0172      0.17121
  2014     2.036      0.11656
  2015    1.9182     0.099836
  2016    2.0025      0.25519
  2017     2.033      0.14581
  2018       NaN          NaN
>>

Other functions can be done similarly...

[Bishwajit Roy]

Hi, When I run this codes, I get the following error:

Error in TrendinYears (line 1) opts=detectImportOptions('barguna.txt');

Second, the mean I receive by doing it in excel is different than you have found. Could it be due to the missing data in my time series? I need to compute the trends (in mm/yr) at each site as well as standard error, to see if the trends are statistically significant or not!! What i will receive by this code probably the mean of all the values/year.

Can you tell me how to solve it? Regards Roy

[dpb]

I shortened the file name locally, open/read the file you have, of course.

Yes, the above is the mean of all finite values for each year.

Use grouping variables to accumulate by other divisions; there is no data available in the dataset provided by which to identify sites so it's not possible to illustrate with that data file alone.

See

doc findgroups
doc splitapply

and examples for doing such calculations. Read the documentation carefully and explore links within it; don't just rely on somebody else to solve all your problems for you.

[dpb]

"used the 12 months running mean to calculate the annual mean" means (so to speak :) ) what, precisely?

What is a 12-month running annual mean--just the yearly average? If so, that's what's returned above; any difference w/ Excel is owing to either data or some other definition of what actually was calculated or, perhaps, a mistake.

[Bishwajit Roy]

Hi, I might have made some mistakes in the explanation of the questions which probably made you confused. And, i used the running mean to remove the seasonal cycle on a monthly basis as you can see from the code provided above. I was not just sure if it's correct. Therefore, wanted to be clear my doubts along with the other questions of calculating the trend and standard error. No problem! You have been very helpful so far! Thanks again!

[dpb]

So you're removing "seasonal" trend by subtracting the mean for each month before computing the annual mean?

That's easily-enough done, but there are probably better ways to handle seasonality unless this is a recognized, standard technique in the specific field.

[Bishwajit Roy]

Hi again, No. At first I am calculating the annual moving average (mean) and then subtracting it from each daily observation to remove the seasonality.

[dpb]

What's the "the annual moving average (mean)", specifically?

[Bishwajit Roy]

removing the annual (or seasonal) cycle by fitting a 12 month running mean. Hope you have understood now.

[dpb]

Not, really, no. :)

What is the definition of that computation as you define it?

Maybe the other way to ask the question is what is the length of the run and over what time period that is subtracted for the given day?

Need a formulation to figure out how to implement it if going to use more "Matlab-like" calculation such as retime

Without resorting to more effort/time than I have to invest, it surely looks to me like you're calculating an average on a month-by-month basis in the above code???

And, I don't see anywhere that that is being subtracted from the individual readings. Where, specifically, is that?

[Bhowmik.U]

what if alternative of retime in matlab earlier versions?

[dpb]

How early?