In this exercise we want to predict daily mean wind speeds in Potsdam. Our model is a simple AR(1) model \[ x_t = a x_{t-1} + v\xi_t, \] where the Gaussian noise xi has unit variance and zero mean. So the parameters to be optimized are 'a' and 'v'.
The data can be found at https://www.ecad.eu/dailydata/customquery.php. Select 'blended', 'Germany', 'Potsdam', and 'wind speed' and download the txt-file. Import the time series in python and plot it. The complete time series has a lot of missing data. We can restrict ourselves to a time window of 10000 points, by selecting data=data[20000:30000]
. Then we can look for the best parameter 'a'.
from scipy.optimize import minimize
def tamsdAR1( DELTA , K , a ):
return 2*K * ( 1 - a**DELTA )
def AR1fit( DELTA, TAMSD ):
def fct( x ):
if(abs(x[1])>0.999):
return 1e45*x[1]**2
return np.mean( np.log(TAMSD / tamsdAR1( DELTA , x[0] , x[1] ))**2 )
x0 = [np.median(TAMSD)/2,0.8]
K,a = minimize( fct , x0 , method='Nelder-Mead' ).x
return K , a
Fit the TAMSD for DELTA smaller than 50, by calling the function as K,a=AR1fit(DELTA[:50],TAMSD[:50])
. Plot the TAMSD and the fit. What is the correlation time of your fitted model?scipy.minimize
function. What is the correlation time of your fitted model?