"Yes, the parameters make sense. The average pressure is about 1013 hPa (the intercept); the slope is 0, so there is no increasing/decreasing trend. Looking at the periodic signals we can see that the daily parameter for theta is 1.6, which is around half of pi, or one-quarter of the daily period; therefore the pressure is highest around 6:00-7:00 AM. Looking at the amplitude parameters we can see that the yearly periodic signal has the highest impact on the pressure.\n",
"Yes, the parameters make sense. The average pressure is about 1013 hPa (the intercept); the slope is 0, so there is no increasing/decreasing trend. Looking at the periodic signals we can see that the daily parameter for theta is 1.6, which is around half of pi, or one-quarter of the daily period. Therefore, the cosine shifts to the left (by a quarter wavelength) and the pressure is lowest around 6 AM and highest around 6 PM (18:00), which would make sense physically for this (simulated) dataset. For the yearly cycles, theta_1 is about - pi, meaning that the cosine shifts to the right (by half a wavelength), meaning the lowest pressure is in Winder (Dec/Jan), and highest pressure in Summer (Jun/Jul). Finally, looking at the amplitude parameters, we can see that the yearly periodic signal has the highest impact on the pressure.\n",
*[CEGM1000 MUDE](http://mude.citg.tudelft.nl/): Week 2.4, Time series analysis. For: December 4, 2024.*
*[CEGM1000 MUDE](http://mude.citg.tudelft.nl/): Week 2.4, Time series analysis. For: December 4, 2024.*
%% Cell type:markdown id: tags:
%% Cell type:markdown id: tags:
In this workshop we will analyze Atmospheric pressure data. The data is a time series of atmospheric pressure measurements over a period of 2 years, with measurements taken every 4 hours. The first column is the time in days, the second column is the atmospheric pressure in hPa.
In this workshop we will analyze Atmospheric pressure data. The data is a time series of atmospheric pressure measurements over a period of 2 years, with measurements taken every 4 hours. The first column is the time in days, the second column is the atmospheric pressure in hPa.
In this workshop we will analyze the time series. We will first fit a functional model to the data in order to stationarize the data. Then we will address the remaining noise and establish a stochastic model by fitting an AR(1) model to the residuals.
In this workshop we will analyze the time series. We will first fit a functional model to the data in order to stationarize the data. Then we will address the remaining noise and establish a stochastic model by fitting an AR(1) model to the residuals.
We begin by loading the data and plotting it.
We begin by loading the data and plotting it.
%% Cell type:markdown id: tags:
%% Cell type:markdown id: tags:
## Part 0.1: Load the data and plot it
## Part 0.1: Load the data and plot it
This week we are using a new library called `statsmodels`. This library is very useful for time series analysis. If you don't have it installed, you can install it by running `!pip install statsmodels` in a code cell. You can also install it by running `pip install statsmodels` in a terminal of your choice (make sure you are in the right environment activated!).
This week we are using a new library called `statsmodels`. This library is very useful for time series analysis. If you don't have it installed, you can install it by running `!pip install statsmodels` in a code cell. You can also install it by running `pip install statsmodels` in a terminal of your choice (make sure you are in the right environment activated!).
%% Cell type:code id: tags:
%% Cell type:code id: tags:
``` python
``` python
import numpy as np
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.pyplot as plt
from statsmodels.graphics.tsaplots import plot_acf
from statsmodels.graphics.tsaplots import plot_acf
from scipy.stats import chi2
from scipy.stats import chi2
from scipy.signal import periodogram
from scipy.signal import periodogram
```
```
%% Cell type:code id: tags:
%% Cell type:code id: tags:
``` python
``` python
data = np.loadtxt('atm_data.txt', delimiter=',')
data = np.loadtxt('atm_data.txt', delimiter=',')
time = data[:, 0]
time = data[:, 0]
data = data[:, 1]
data = data[:, 1]
# dt = YOUR_CODE_HERE # Time step
# dt = YOUR_CODE_HERE # Time step
# fs = YOUR_CODE_HERE # Sampling frequency
# fs = YOUR_CODE_HERE # Sampling frequency
# SOLUTION
# SOLUTION
dt = time[1] - time[0]
dt = time[1] - time[0]
fs = 1 / dt
fs = 1 / dt
plt.figure(figsize=(10, 3))
plt.figure(figsize=(10, 3))
plt.plot(time, data)
plt.plot(time, data)
plt.xlabel('Time [days]')
plt.xlabel('Time [days]')
plt.ylabel('Atmospheric Pressure [hPa]')
plt.ylabel('Atmospheric Pressure [hPa]')
plt.title('2 year of atmospheric pressure data')
plt.title('2 year of atmospheric pressure data')
plt.grid(True)
plt.grid(True)
```
```
%% Output
%% Output
%% Cell type:markdown id: tags:
%% Cell type:markdown id: tags:
## Part 1: Finding frequencies of periodic patterns
## Part 1: Finding frequencies of periodic patterns
We clearly see that the data contains a seasonal pattern. We will start by fitting a functional model to the data in order to make the data stationary. To find the frequency of the seasonal pattern we will use the power spectrum of the data.
We clearly see that the data contains a seasonal pattern. We will start by fitting a functional model to the data in order to make the data stationary. To find the frequency of the seasonal pattern we will use the power spectrum of the data.
The time series may be influenced by multiple periodic signals at different frequencies; therefore we will apply a method to find the dominant frequency, and then iteratively remove the corresponding signal from the data.
The time series may be influenced by multiple periodic signals at different frequencies; therefore we will apply a method to find the dominant frequency, and then iteratively remove the corresponding signal from the data.
We will create a function `find_frequency` that will take the data and an $\mathrm{A}$-matrix as input and return the frequency of the most dominant periodic signal. Note that we cannot have an empty A-matrix when running the function for the first time. Therefore, we first define an A matrix for a linear model that accounts for an intercept and slope, with which `find_frequency` will then de-trend the time series.
We will create a function `find_frequency` that will take the data and an $\mathrm{A}$-matrix as input and return the frequency of the most dominant periodic signal. Note that we cannot have an empty A-matrix when running the function for the first time. Therefore, we first define an A matrix for a linear model that accounts for an intercept and slope, with which `find_frequency` will then de-trend the time series.
The function will look like this:
The function will look like this:
1. Using the $\mathrm{A}$-matrix we will fit a model to the data using the least squares method.
1. Using the $\mathrm{A}$-matrix we will fit a model to the data using the least squares method.
2. We will calculate the PSD of the residuals of the model.
2. We will calculate the PSD of the residuals of the model.
3. We will find the frequency of the most dominant periodic signal in the PSD.
3. We will find the frequency of the most dominant periodic signal in the PSD.
4. Optional: we can plot the power spectrum of the data with the highest peak.
4. Optional: we can plot the power spectrum of the data with the highest peak.
Read the code in the next Python cell study the contents until you understand what the functions are doing. Then complete the missing parts of the code.
Read the code in the next Python cell study the contents until you understand what the functions are doing. Then complete the missing parts of the code.
We will write functions that return multiple parameters, but if we are not interested in defining a particular parameter we can use `_` to ignore it. For example, if we write <code>_, b = function(a)</code>, we only define the second output of the function; the first output is not stored in a variable.</p></div>
We will write functions that return multiple parameters, but if we are not interested in defining a particular parameter we can use `_` to ignore it. For example, if we write <code>_, b = function(a)</code>, we only define the second output of the function; the first output is not stored in a variable.</p></div>
%% Cell type:code id: tags:
%% Cell type:code id: tags:
``` python
``` python
def fit_model(data, time, A, plot=False):
def fit_model(data, time, A, plot=False):
'''
'''
Function to find the least squares solution of the data
Function to find the least squares solution of the data
data: input data
data: input data
time: time vector
time: time vector
A: A-matrix to fit the data
A: A-matrix to fit the data
plot: boolean to plot the results or not
plot: boolean to plot the results or not
'''
'''
# x_hat = YOUR_CODE_HERE # least squares solution
# x_hat = YOUR_CODE_HERE # least squares solution
# y_hat = YOUR_CODE_HERE # model prediction
# y_hat = YOUR_CODE_HERE # model prediction
# e_hat = YOUR_CODE_HERE # residuals
# e_hat = YOUR_CODE_HERE # residuals
# SOLUTION
# SOLUTION
x_hat = np.linalg.solve(A.T @ A, A.T @ data)
x_hat = np.linalg.solve(A.T @ A, A.T @ data)
y_hat = A @ x_hat
y_hat = A @ x_hat
e_hat = data - y_hat
e_hat = data - y_hat
# END SOLUTION BLOCK
# END SOLUTION BLOCK
if plot:
if plot:
plt.figure(figsize=(10, 5))
plt.figure(figsize=(10, 5))
plt.subplot(211)
plt.subplot(211)
plt.plot(time, data, label='Data')
plt.plot(time, data, label='Data')
plt.plot(time, y_hat, label='Estimated data')
plt.plot(time, y_hat, label='Estimated data')
plt.xlabel('Time [days]')
plt.xlabel('Time [days]')
plt.ylabel('Atmospheric Pressure [hPa]')
plt.ylabel('Atmospheric Pressure [hPa]')
plt.title('Data vs Estimated data')
plt.title('Data vs Estimated data')
plt.grid(True)
plt.grid(True)
plt.legend()
plt.legend()
plt.subplot(212)
plt.subplot(212)
plt.plot(time, e_hat, label='Residuals')
plt.plot(time, e_hat, label='Residuals')
plt.xlabel('Time [days]')
plt.xlabel('Time [days]')
plt.ylabel('Atmospheric Pressure [hPa]')
plt.ylabel('Atmospheric Pressure [hPa]')
plt.title('Residuals')
plt.title('Residuals')
plt.grid(True)
plt.grid(True)
plt.legend()
plt.legend()
plt.tight_layout()
plt.tight_layout()
return x_hat, y_hat, e_hat
return x_hat, y_hat, e_hat
def find_frequency(data, time, A, fs, plot=True):
def find_frequency(data, time, A, fs, plot=True):
'''
'''
Function to find the dominant frequency of the signal
Function to find the dominant frequency of the signal
data: input data
data: input data
time: time vector
time: time vector
A: A-matrix to detrend the data (prior to spectral analysis)
A: A-matrix to detrend the data (prior to spectral analysis)
How can we use the <code>find_frequency</code> function to find the frequency of all periodic pattern in the data? In other words, how can we iteratively detect the frequencies in the data?
How can we use the <code>find_frequency</code> function to find the frequency of all periodic pattern in the data? In other words, how can we iteratively detect the frequencies in the data?
Write your answer in a bulleted list that describes the procedure in the Markdown cell below.
Write your answer in a bulleted list that describes the procedure in the Markdown cell below.
<em>Hint: we should try to remove frequencies from the data once we have detected them.</em>
<em>Hint: we should try to remove frequencies from the data once we have detected them.</em>
Now we will run `find_frequency` repeatedly to find the frequencies in the data. For the first iteration we will use an $\mathrm{A}$-matrix that will remove a linear trend from the data. When we have found a frequency we will also remove the corresponding signal from the data in the next iteration.
Now we will run `find_frequency` repeatedly to find the frequencies in the data. For the first iteration we will use an $\mathrm{A}$-matrix that will remove a linear trend from the data. When we have found a frequency we will also remove the corresponding signal from the data in the next iteration.
Look back at the plot of the original time series and confirm whether or not you can see the presence of the signal with the frequency found in the previous task. Write your answer in one sentence in the Markdown cell below.
Look back at the plot of the original time series and confirm whether or not you can see the presence of the signal with the frequency found in the previous task. Write your answer in one sentence in the Markdown cell below.
Find a second frequency by repeating the process. Before using <code>find_frequency</code> again, the linear trend and first periodic component should be accounted for.
Find a second frequency by repeating the process. Before using <code>find_frequency</code> again, the linear trend and first periodic component should be accounted for.
Repeat this step until you have removed all dominant signals (you will have to determine when to stop).
Repeat this step until you have removed all dominant signals (you will have to determine when to stop).
Print the frequencies of all periodic signals found in the data.
Print the frequencies of all periodic signals found in the data.
</p>
</p>
</div>
</div>
%% Cell type:code id: tags:
%% Cell type:code id: tags:
``` python
``` python
# YOUR_CODE_HERE # may be more than one line or more than one cell
# YOUR_CODE_HERE # may be more than one line or more than one cell
Describe the number of steps taken in the previous task and explain in your own words how the dominant frequencies were determined (print them, if you did not already do so).
Describe the number of steps taken in the previous task and explain in your own words how the dominant frequencies were determined (print them, if you did not already do so).
How did you decide when to stop?
How did you decide when to stop?
Is the final (detrended) time series (residuals) stationary? Explain.
Is the final (detrended) time series (residuals) stationary? Explain.
- First run: we fit a linear model to the data and find the dominant frequency of 1 cycle per year. PSD of the residuals are not stationary, so we need to do this again.
- First run: we fit a linear model to the data and find the dominant frequency of 1 cycle per year. PSD of the residuals are not stationary, so we need to do this again.
- Second run: we fit a model with a linear trend and a yearly cycle. We find a dominant frequency of 1 cycle per day.
- Second run: we fit a model with a linear trend and a yearly cycle. We find a dominant frequency of 1 cycle per day.
- Third run: the highest PSD is a signal with frequency 0.07 cycle/day, which is about 3 orders of magnitude less than the previous run (comparing the order of magnitude on the y-axis between the previous two plots). Therefore, the third frequency is not necessary.
- Third run: the highest PSD is a signal with frequency 0.07 cycle/day, which is about 3 orders of magnitude less than the previous run (comparing the order of magnitude on the y-axis between the previous two plots). Therefore, the third frequency is not necessary.
- The previous point illustrates that the detrended series is indeed stationary.
- The previous point illustrates that the detrended series is indeed stationary.
</p>
</p>
</div>
</div>
%% Cell type:markdown id: tags:
%% Cell type:markdown id: tags:
## Part 2: Fitting the functional model
## Part 2: Fitting the functional model
In the next cell we will fit the model to generate stationary residuals. Above we have used for each dominant frequency $f_i$ ($i=1,2$) the model:
In the next cell we will fit the model to generate stationary residuals. Above we have used for each dominant frequency $f_i$ ($i=1,2$) the model:
$$a_i \cos(2\pi f_i t) + b_i \sin(2\pi f_i t)$$
$$a_i \cos(2\pi f_i t) + b_i \sin(2\pi f_i t)$$
However, to report the periodic signals we would like to have the amplitude, phase shift and the frequency of those signals, which can be recovered from:
However, to report the periodic signals we would like to have the amplitude, phase shift and the frequency of those signals, which can be recovered from:
$$A_i \cos(2\pi f_i t + \theta_i)$$
$$A_i \cos(2\pi f_i t + \theta_i)$$
Where the amplitude $A_i = \sqrt{a_i^2 + b_i^2}$ and $\theta_i = \arctan(-b_i/a_i)$
Where the amplitude $A_i = \sqrt{a_i^2 + b_i^2}$ and $\theta_i = \arctan(-b_i/a_i)$
Note: in Section 4.1 book this was shown where the angular frequency $\omega = 2\pi f$ was used.
Note: in Section 4.1 book this was shown where the angular frequency $\omega = 2\pi f$ was used.
If you want to compute the $\theta_i$ you can use the <code>np.arctan2</code> function. This function is a version of the arctan function that returns the angle in the correct quadrant. Using the <code>np.arctan</code> function will not give you the correct angle!</p></div>
If you want to compute the $\theta_i$ you can use the <code>np.arctan2</code> function. This function is a version of the arctan function that returns the angle in the correct quadrant. Using the <code>np.arctan</code> function will not give you the correct angle!</p></div>
%% Cell type:code id: tags:
%% Cell type:code id: tags:
``` python
``` python
def rewrite_seasonal_comp(ak, bk):
def rewrite_seasonal_comp(ak, bk):
'''
'''
Function to rewrite the seasonal component in terms of sin and cos
Function to rewrite the seasonal component in terms of sin and cos
<div style="background-color:#facb8e; color: black; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px; width: 95%"> <p>Task 4.2 was moved here, since it made more sense to interpret the parameters earlier.</p></div>
<div style="background-color:#facb8e; color: black; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px; width: 95%"> <p>Task 4.2 was moved here, since it made more sense to interpret the parameters earlier.</p></div>
Yes, the parameters make sense. The average pressure is about 1013 hPa (the intercept); the slope is 0, so there is no increasing/decreasing trend. Looking at the periodic signals we can see that the daily parameter for theta is 1.6, which is around half of pi, or one-quarter of the daily period; therefore the pressure is highest around 6:00-7:00 AM. Looking at the amplitude parameters we can see that the yearly periodic signal has the highest impact on the pressure.
Yes, the parameters make sense. The average pressure is about 1013 hPa (the intercept); the slope is 0, so there is no increasing/decreasing trend. Looking at the periodic signals we can see that the daily parameter for theta is 1.6, which is around half of pi, or one-quarter of the daily period. Therefore, the cosine shifts to the left (by a quarter wavelength) and the pressure is lowest around 6 AM and highest around 6 PM (18:00), which would make sense physically for this (simulated) dataset. For the yearly cycles, theta_1 is about - pi, meaning that the cosine shifts to the right (by half a wavelength), meaning the lowest pressure is in Winder (Dec/Jan), and highest pressure in Summer (Jun/Jul). Finally, looking at the amplitude parameters, we can see that the yearly periodic signal has the highest impact on the pressure.
</p>
</p>
</div>
</div>
%% Cell type:markdown id: tags:
%% Cell type:markdown id: tags:
## Part 3: Analyzing the residuals
## Part 3: Analyzing the residuals
Now that we have our residuals we can fit an AR model to the residuals. We will start by plotting the ACF of the residuals. We will then fit an AR model to the residuals and report the parameters of the AR model..
Now that we have our residuals we can fit an AR model to the residuals. We will start by plotting the ACF of the residuals. We will then fit an AR model to the residuals and report the parameters of the AR model..
The data is autocreated since lags are significant correlated (up to lag 3)
The data is autocreated since lags are significant correlated (up to lag 3)
</p>
</p>
</div>
</div>
%% Cell type:markdown id: tags:
%% Cell type:markdown id: tags:
### Fit an AR(1) model to the residuals
### Fit an AR(1) model to the residuals
First we write a function `AR1` that will take the stationary residuals as input and return the parameters of the AR1 model. Then we will fit this model to the residuals and report the parameters.
First we write a function `AR1` that will take the stationary residuals as input and return the parameters of the AR1 model. Then we will fit this model to the residuals and report the parameters.
# plt.xlim([0, 100]) # uncomment this line to zoom in, for better visualization
# plt.xlim([0, 100]) # uncomment this line to zoom in, for better visualization
plt.grid(True)
plt.grid(True)
plt.legend()
plt.legend()
print(f'Estimated Parameters:')
print(f'Estimated Parameters:')
print(f'phi = {x_hat[0]:.4f}')
print(f'phi = {x_hat[0]:.4f}')
return x_hat, e_hat
return x_hat, e_hat
```
```
%% Cell type:markdown id: tags:
%% Cell type:markdown id: tags:
### Fit an AR1 model to the residuals
### Fit an AR1 model to the residuals
Fit an AR(1) model to the residuals of the functional model (found in Part 2), then report the parameters. We will need to check if the AR(1) model is a good fit to the residuals, although there are many ways to check whether this is the case. We can use a testing procedure to compare different AR models, an example of this can be found in the book.
Fit an AR(1) model to the residuals of the functional model (found in Part 2), then report the parameters. We will need to check if the AR(1) model is a good fit to the residuals, although there are many ways to check whether this is the case. We can use a testing procedure to compare different AR models, an example of this can be found in the book.
Today is your lucky day: we will use a more simple method. If the AR(1) model is a good fit to the residuals, the residuals of the AR(1) model (output by the `AR1` function) should be white noise. We will plot the ACF of the residuals of the AR(1) model to check if this is the case.
Today is your lucky day: we will use a more simple method. If the AR(1) model is a good fit to the residuals, the residuals of the AR(1) model (output by the `AR1` function) should be white noise. We will plot the ACF of the residuals of the AR(1) model to check if this is the case.
No autocorrelation in the residuals of the AR(1) model, so the AR(1) model is a good fit to the residuals. If the AR(1) model is not a good fit to the residuals we could try a higher order AR model.
No autocorrelation in the residuals of the AR(1) model, so the AR(1) model is a good fit to the residuals. If the AR(1) model is not a good fit to the residuals we could try a higher order AR model.
</p>
</p>
</div>
</div>
%% Cell type:markdown id: tags:
%% Cell type:markdown id: tags:
## Part 4: Report the results
## Part 4: Report the results
_This part is optional because it is not part of the book, and is not part of the exam material. However, it is a nice exercise to see how the uncertainty in our model parameters can be computed for situations where the observations are **not** independent. This is different from Q1 when we studied observation theory models with an assumption of independent observations. Maybe we will also (optionally) use this on Friday..._
_This part is optional because it is not part of the book, and is not part of the exam material. However, it is a nice exercise to see how the uncertainty in our model parameters can be computed for situations where the observations are **not** independent. This is different from Q1 when we studied observation theory models with an assumption of independent observations. Maybe we will also (optionally) use this on Friday..._
Now that we have found the periodic signals in the data and fitted an AR model to the residuals, we can report the results. By combining including the AR (noise) process, we get residuals that are white noise. When the model hat white noise residuals, we can also report the confidence intervals of the model.
Now that we have found the periodic signals in the data and fitted an AR model to the residuals, we can report the results. By combining including the AR (noise) process, we get residuals that are white noise. When the model hat white noise residuals, we can also report the confidence intervals of the model.
We will use the unbiased estimate of the variance of the residuals to calculate the confidence intervals. The unbiased estimate of the variance is given by:
We will use the unbiased estimate of the variance of the residuals to calculate the confidence intervals. The unbiased estimate of the variance is given by:
<div style="background-color:#facb8e; color: black; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px; width: 95%"> <p>Task 4.2 was moved above in the notebook, since it made more sense to interpret the parameters earlier.</p></div>
<div style="background-color:#facb8e; color: black; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px; width: 95%"> <p>Task 4.2 was moved above in the notebook, since it made more sense to interpret the parameters earlier.</p></div>
Yes, the parameters make sense. The average pressure is about 1013 hPa (the intercept); the slope is 0, so there is no increasing/decreasing trend. Looking at the periodic signals we can see that the daily parameter for theta is 1.6, which is around half of pi, or one-quarter of the daily period; therefore the pressure is highest around 6:00-7:00 AM. Looking at the amplitude parameters we can see that the yearly periodic signal has the highest impact on the pressure.
Yes, the parameters make sense. The average pressure is about 1013 hPa (the intercept); the slope is 0, so there is no increasing/decreasing trend. Looking at the periodic signals we can see that the daily parameter for theta is 1.6, which is around half of pi, or one-quarter of the daily period. Therefore, the cosine shifts to the left (by a quarter wavelength) and the pressure is lowest around 6 AM and highest around 6 PM (18:00), which would make sense physically for this (simulated) dataset. For the yearly cycles, theta_1 is about - pi, meaning that the cosine shifts to the right (by half a wavelength), meaning the lowest pressure is in Winder (Dec/Jan), and highest pressure in Summer (Jun/Jul). Finally, looking at the amplitude parameters, we can see that the yearly periodic signal has the highest impact on the pressure.
