pymcssa

  5class AR1estimator:
  6    def __init__(self, data):
  7        self.data = np.asarray(data)
  8        self.n = len(data)
  9    def c_hats(self,l):
 10        """Estimate lag-l sample covariance (biased by n-l).
 11
 12        Computes c_hat(l) = sum_{i=0}^{n-l-1} (x_i - mean)*(x_{i+l} - mean) / (n-l)
 13
 14        Args:
 15            l (int): Lag value (non-negative).
 16
 17        Returns:
 18            float: Estimated covariance at lag l.
 19        """
 20        ch=0
 21        for i in range(len(self.data)-l):
 22            ch+=((self.data[i]-np.mean(self.data))*(self.data[i+l]-np.mean(self.data)))
 23        ch=ch/(self.n-l)
 24        return ch
 25    def mu2(self,gam):
 26        """Compute mu2 function used in AR(1) parameter estimation.
 27
 28        Args:
 29            gam (float): AR(1) coefficient.
 30
 31        Returns:
 32            float: Value of mu2(gam) used in Newton-Raphson iterations.
 33        """
 34        mu2=0        
 35        for i in range(self.n-1):
 36            mu2+=2*(self.n-(i+1))*(gam**(i+1))
 37        mu2=(1/self.n)+((1/(self.n*self.n))*mu2)
 38        return mu2
 39    def tmumudash(self,gamma):
 40        """Compute derivative-like term tmumudash used in nr_dash.
 41
 42        Args:
 43            gamma (float): AR(1) coefficient.
 44
 45        Returns:
 46            float: Value used to construct derivative in nr_dash.
 47        """
 48        mud=0
 49        for i in range(self.n-1):
 50            mud+=(i+1)*(self.n-(i+1))*(gamma**i)
 51        mud=2*mud/(self.n*self.n)
 52        return mud
 53    def nr_func(self,gamma):
 54        """Newton-Raphson objective function for estimating gamma.
 55
 56        The equation solved is (gamma - mu2(gamma)) / (1 - mu2(gamma)) = c1/c0.
 57
 58        Args:
 59            gamma (float): Current gamma estimate.
 60
 61        Returns:
 62            float: Function value.
 63
 64        Raises:
 65            ZeroDivisionError: If estimated c0 (lag 0) is zero.
 66        """
 67        c1=self.c_hats(1)
 68        c0=self.c_hats(0)
 69        if c0 == 0:
 70            raise ZeroDivisionError("c0 is zero in nr_func")
 71        f=((gamma-self.mu2(gamma))/(1-self.mu2(gamma)))-(c1/c0)
 72        return f
 73    def nr_dash(self,gamma):
 74        """Approximate derivative of the NR objective for gamma update.
 75
 76        Args:
 77            gamma (float): Current gamma estimate.
 78
 79        Returns:
 80            float: Derivative approximation used in NR step.
 81
 82        Raises:
 83            ZeroDivisionError: If estimated c0 (lag 0) is zero.
 84        """
 85        c1=self.c_hats(1)
 86        c0=self.c_hats(0)
 87        if c0 == 0:
 88            raise ZeroDivisionError("c0 is zero in nr_dash")
 89        denom=(1-self.mu2(gamma))
 90        if denom == 0:
 91            denom=1e-10
 92        fd=1+(self.tmumudash(gamma)*(self.nr_func(gamma)+(c1/c0)-1))
 93        fd=fd/denom
 94        return fd
 95    def gambar(self,max_iter=1000):
 96        """Estimate AR(1) coefficient gamma via Newton-Raphson.
 97
 98        Starts from c1/c0 and iterates until relative error < 1e-6 or max_iter reached.
 99
100        Args:
101            max_iter (int): Maximum NR iterations.
102
103        Returns:
104            float: Estimated gamma value (also stored as self.gam).
105        """
106        c1=self.c_hats(1)
107        c0=self.c_hats(0)
108        ini=c1/c0
109        err=100
110        iter_count=0
111        while err>1e-6 and iter_count<max_iter:
112            ini_prev=ini
113            nr_f=self.nr_func(ini)
114            nr_d=self.nr_dash(ini)
115            if nr_d == 0:
116                break
117            ini=ini-(nr_f/nr_d)
118            err=abs((ini-ini_prev)/(ini_prev if ini_prev!= 0 else 1))
119            iter_count += 1
120        self.gam=ini
121        return ini
122    def alph(self,max_iter=1000):
123        """Estimate AR(1) innovation standard deviation alpha.
124
125        Uses the estimated gamma and mu2 to compute alpha = sqrt(c0bar * (1 - gamma)),
126        where c0bar = c0 / (1 - mu2(gamma)). Small positive floor applied to avoid
127        negative or zero divisors.
128
129        Args:
130            max_iter (int): Maximum NR iterations.
131
132        Returns:
133            float: Estimated alpha (also stored as self.alpha).
134        """
135        gam=self.gambar(max_iter=max_iter)
136        c0=self.c_hats(0)
137        mu2_val=self.mu2(gam)
138        denom=(1-mu2_val)
139        if denom <= 0:
140            denom=1e-10
141        c0bar=c0/denom
142        sqrt_term=c0bar*(1-gam)
143        sqrt_term=sqrt_term if sqrt_term > 0 else 1e-10  # prevent sqrt of negative
144        alpha=ma.sqrt(sqrt_term)
145        self.alpha=alpha
146        return alpha
147
148    def ar1_model(self,alpha,gamma):
149        """Generate a single AR(1) surrogate time series of length n.
150
151        Model: r[t] = gamma * r[t-1] + alpha * w[t] with w ~ N(0,1) and r[0]=0.
152
153        Args:
154            alpha (float): Innovation std.
155            gamma (float): AR(1) coefficient.
156
157        Returns:
158            ndarray: 1-D surrogate time series of length n.
159        """
160        rn=np.zeros(self.n)
161        wn=np.random.normal(size=self.n)
162        for i in range(1,self.n):
163            rn[i]=gamma*(rn[i-1])+alpha*wn[i]
164        return rn

  6class SSA:
  7    """
  8    Basic Single-Channel Singular Spectrum Analysis (SSA).
  9
 10    Performs Singular Spectrum Analysis (SSA) decomposition and reconstruction
 11    of a 1D time series using Singular Value Decomposition (SVD) of the
 12    trajectory (Hankel) matrix.
 13
 14    This method decomposes the input signal into interpretable components such
 15    as trends, oscillations, and noise, based on temporal correlations.
 16
 17    Attributes
 18    ----------
 19    data : ndarray
 20        Input 1D time series.
 21    m : int
 22        Embedding window length.
 23    k_vals : int
 24        Number of leading components to reconstruct.
 25
 26    Notes
 27    -----
 28    SSA is a non-parametric spectral decomposition technique that extracts
 29    dominant temporal patterns without assuming any underlying model.
 30    """
 31    def __init__(self,data,m,k_vals):
 32        """
 33        Initialize the SSA object.
 34
 35        Parameters
 36        ----------
 37        data : array-like
 38            Input 1D time series.
 39        m : int
 40            Embedding (window) length.
 41        k_vals : int
 42            Number of reconstructed components to retain.
 43        """
 44        self.data=np.asarray(data)
 45        self.m=m
 46        self.k_vals=k_vals
 47    def ssa(self):
 48        """
 49        Perform SSA decomposition and reconstruction.
 50
 51        Constructs the trajectory matrix using `Tmat`, performs SVD to obtain
 52        eigenvectors (EOFs) and singular values, computes principal components,
 53        and reconstructs the leading `k_vals` components.
 54
 55        Returns
 56        -------
 57        results : dict
 58            Dictionary containing:
 59            
 60            - **eofs** : ndarray, shape (m, m)
 61                EOFs (Empirical Orthogonal Functions) or spatial patterns.
 62            - **eigen_values** : ndarray, shape (m,)
 63                Eigenvalues (squared singular values) representing variance.
 64            - **percent_explained** : ndarray, shape (m,)
 65                Percentage of variance explained by each component.
 66            - **pcs** : ndarray, shape (n - m + 1, m)
 67                Principal components (temporal coefficients).
 68            - **rcs** : ndarray, shape (n, k_vals)
 69                Reconstructed time series components using first `k_vals` modes.
 70
 71        Raises
 72        ------
 73        Exception
 74            If the input data is multi-dimensional (use MSSA or PCA instead).
 75        """
 76
 77        n=np.shape(self.data)[0]
 78        if len(np.shape(self.data))>1:
 79            raise Exception("multi dimensional data use mssa or pca. can't make an ssa here")
 80        elif len(np.shape(self.data))==1:
 81            nn=n-self.m+1
 82            T=Tmat(self.data,self.m)
 83            T=T*(1/ma.sqrt(nn))
 84            T=np.array(T)
 85            eofs,s,rhot=np.linalg.svd(T)
 86            pcs=np.zeros((nn,self.m))
 87            for t in range(nn):
 88                for k in range(self.m):
 89                    for j in range(self.m):
 90                        pcs[t,k]+=(self.data[t+j]*eofs[j,k])
 91            p_vars=s*s
 92            exps=(p_vars/np.sum(p_vars))*100
 93            params=[0.0,0.0,0.0]
 94            rcs=np.zeros((n,self.k_vals))
 95            for k in range(self.k_vals):
 96                for t in range(n):
 97                    if t>=0 and t<=(self.m-2):
 98                        params[0]=(t+1)
 99                        params[1]=0
100                        params[2]=t+1
101                    if t>=self.m-1 and t<=nn-1:
102                        params[0]=self.m
103                        params[1]=0
104                        params[2]=self.m
105                    if t>=nn and t<=n-1:
106                        params[0]=(n-t)
107                        params[1]=t-n+self.m
108                        params[2]=self.m
109                    for j in range(params[1],params[2]):
110                        rcs[t,k]=rcs[t,k]+((1/params[0])*(pcs[t-j,k]*eofs[j,k]))
111            results={'eofs':eofs,
112                     'eigen_values':p_vars,
113                     'percent_explained':exps,
114                     'pcs':pcs,
115                     'rcs':rcs}
116            return results

 19class MCSSA(AR1estimator):
 20    """
 21    Monte Carlo Singular Spectrum Analysis (MCSSA).
 22
 23    Extends AR1estimator to provide MCSSA workflows that generate AR(1)
 24    surrogate ensembles, compute trajectory/covariance matrices, and produce
 25    percentile-based confidence bounds for SSA eigenvalues.
 26
 27    Usage:
 28        mc = MCSSA(data, m)
 29        results = mc.mcssa_basic(up_perc=97.5, down_perc=2.5, ns=1000)
 30
 31    Attributes:
 32        data (ndarray): 1-D input time series as numpy array.
 33        m (int): Window length for SSA embedding.
 34        n (int): Length of the input time series.
 35        nd (int): Number of rows in trajectory matrix (n - m + 1).
 36        gamma (float): Estimated AR(1) coefficient (set by inherited methods).
 37        alpha (float): Estimated AR(1) innovation std (set by inherited methods).
 38        surrs (ndarray): Last generated surrogate ensemble (ns x n).
 39    """
 40    def __init__(self,data,m):
 41        """
 42        Initialize MCSSA instance.
 43
 44        Args:
 45            data (array-like): 1-D input time series.
 46            m (int): Window length (embedding dimension). Must satisfy 1 <= m <= len(data).
 47
 48        Raises:
 49            ValueError: If m is larger than the length of the data.
 50        """
 51        super().__init__(data)
 52        self.data=np.asarray(data)
 53        self.m=m
 54        self.n=len(data)
 55        self.nd=self.n-self.m+1
 56        if self.nd<=0:
 57            raise ValueError("m must be<=length of data")
 58    def mcssa_basic(self,up_perc,down_perc,ns=100000,max_iter=1000,return_surrogates=False):
 59        """
 60        Monte-Carlo SSA using the data EOF basis.
 61
 62        Generates `ns` AR(1) surrogates using the AR(1) parameters estimated
 63        from the data, projects each surrogate covariance matrix into the data's
 64        EOF basis and computes upper/lower percentile confidence bounds for the
 65        projected eigenvalues.
 66
 67        Args:
 68            up_perc (float): Upper percentile (e.g. 97.5).
 69            down_perc (float): Lower percentile (e.g. 2.5).
 70            ns (int, optional): Number of surrogate realizations. Default 100000.
 71            max_iter (int, optional): Max iterations for AR(1) estimation routines.
 72            return_surrogates (bool, optional): If True, returned dict includes
 73                the surrogate ensemble under "surrogates".
 74
 75        Returns:
 76            dict: {
 77                "data_eigenvalues": ndarray (sorted descending),
 78                "upper_confidence": list of upper percentile values,
 79                "lower_confidence": list of lower percentile values,
 80                "data_eigenvectors": ndarray (columns are EOFs),
 81                "gamma": float,
 82                "alpha": float,
 83                "spreads": ndarray (ns x m) of projected surrogate variances,
 84                ("surrogates": ndarray) optional
 85            }
 86
 87        Raises:
 88            ValueError: If ns < 1, percentiles out of [0,100], or down_perc >= up_perc.
 89        """
 90        if ns < 1:
 91            raise ValueError("ns must be >= 1")
 92        if not (0 <= down_perc <= 100 and 0 <= up_perc <= 100):
 93            raise ValueError("Percentiles must be between 0 and 100")
 94        if down_perc >= up_perc:
 95            raise ValueError("down_perc must be less than up_perc")
 96        T_dat=Tmat(self.data,self.m)
 97        c=(1/self.nd)*np.transpose(T_dat)@T_dat
 98        evs,eofs=np.linalg.eigh(c)
 99        idx=np.argsort(evs)[::-1]
100        evs=evs[idx]
101        eofs=eofs[:, idx]
102        self.gamma=round(self.gambar(max_iter=max_iter),8)
103        self.alpha=round(self.alph(max_iter=max_iter),8)
104        surrs=[]
105        for i in range(ns):
106            surrs.append(self.ar1_model(self.alpha,self.gamma))
107        surrs=np.vstack(surrs)
108        surrs=surrs-np.mean(surrs,axis=1,keepdims=True)
109        self.surrs=surrs
110        spreads=[]
111        for i in range(np.shape(surrs)[0]):
112            Tm=Tmat(surrs[i,:],self.m)
113            c_s=(1/self.nd)*np.transpose(Tm)@Tm
114            c_s_eval=np.transpose(eofs)@c_s@eofs
115            c_s_eval=np.diag(c_s_eval)
116            spreads.append(c_s_eval)
117        spreads=np.vstack(spreads)
118        ups=[]
119        downs=[]
120        for i in range(np.shape(spreads)[1]):
121            ups.append(np.percentile(spreads[:,i],up_perc))
122            downs.append(np.percentile(spreads[:,i],down_perc))
123        result= {"data_eigenvalues": evs,
124                  "upper_confidence": ups,
125                  "lower_confidence": downs,
126                  "data_eigenvectors": eofs,
127                  "gamma": self.gamma,
128                  "alpha": self.alpha,
129                  "spreads": spreads,
130                  "surrogates": self.surrs,
131                  }
132        if not return_surrogates:
133            result.pop("surrogates")
134        return result
135    def mcssa_ensemble(self,up_perc,down_perc,ns=100000,max_iter=1000,return_surrogates=False):
136        """
137        Ensemble MCSSA using the mean surrogate covariance eigenbasis.
138
139        Steps:
140        - generate ns AR(1) surrogates,
141        - compute each surrogate trajectory covariance and form their mean,
142        - compute the eigenbasis of the mean covariance,
143        - project surrogate covariances and the data covariance onto that basis,
144        - compute percentile confidence bounds from projected surrogate spreads.
145
146        Args:
147            up_perc (float): Upper percentile for confidence.
148            down_perc (float): Lower percentile for confidence.
149            ns (int, optional): Number of surrogate realizations.
150            max_iter (int, optional): Max iterations for AR(1) estimation.
151            return_surrogates (bool, optional): If True include surrogates in result.
152
153        Returns:
154            dict: {
155                "data_eigenvalues": ndarray (projected onto mean surrogate basis),
156                "upper_confidence": list,
157                "lower_confidence": list,
158                "mean_surrogate_eigenvalues": ndarray,
159                "mean_surrogate_eigenvectors": ndarray,
160                "spreads": ndarray (ns x m),
161                "alpha": float,
162                "gamma": float,
163                ("surrogates": ndarray) optional
164            }
165
166        Raises:
167            ValueError: If ns < 1 or invalid percentile arguments.
168        """
169        if ns < 1:
170            raise ValueError("ns must be >= 1")
171        if not (0 <= down_perc <= 100 and 0 <= up_perc <= 100):
172            raise ValueError("Percentiles must be between 0 and 100")
173        if down_perc >= up_perc:
174            raise ValueError("down_perc must be less than up_perc")
175        self.gamma=round(self.gambar(max_iter=max_iter),8)
176        self.alpha=round(self.alph(max_iter=max_iter),8)
177        surrs=[]
178        for i in range(ns):
179            surrs.append(self.ar1_model(self.alpha,self.gamma))
180        surrs=np.vstack(surrs)
181        surrs=surrs-np.mean(surrs,axis=1,keepdims=True)
182        self.surrs=surrs
183        cs=[]
184        for i in range(np.shape(surrs)[0]):
185            Ts=Tmat(surrs[i,:],self.m)
186            cs.append((1/self.nd)*np.transpose(Ts)@Ts)
187        c_m=np.mean(np.stack(cs, axis=0), axis=0)
188        c_m_eval, c_m_evs=np.linalg.eigh(c_m)
189        self.c_mean_eval=c_m_eval
190        self.c_mean_evs=c_m_evs
191        idx=np.argsort(c_m_eval)[::-1]
192        c_m_eval=c_m_eval[idx]
193        c_m_evs=c_m_evs[:, idx]
194        spreads=[]
195        for i in range(np.shape(surrs)[0]):
196            Traj=Tmat(surrs[i],self.m)
197            c_s=(1/self.nd)*np.transpose(Traj)@Traj
198            c_s_eval=c_m_evs.T@c_s@c_m_evs
199            c_s_eval=np.diag(c_s_eval)
200            spreads.append(c_s_eval)
201        spreads=np.vstack(spreads)
202        dats=self.data-np.mean(self.data)
203        T_dat=Tmat(dats,self.m)
204        c=(1/self.nd)*np.transpose(T_dat)@T_dat
205        c_eval=np.transpose(c_m_evs)@c@c_m_evs
206        c_eval=np.diag(c_eval)
207        ups=[]
208        downs=[]
209        for i in range(np.shape(spreads)[1]):
210            ups.append(np.percentile(spreads[:,i],up_perc))
211            downs.append(np.percentile(spreads[:,i],down_perc))
212        results={"data_eigenvalues": c_eval,
213                   "upper_confidence": ups,
214                   "lower_confidence": downs,
215                   "mean_surrogate_eigenvalues": c_m_eval,
216                   "mean_surrogate_eigenvectors": c_m_evs,
217                   "surrogates": surrs,
218                   "spreads": spreads,
219                   "alpha": self.alpha,
220                   "gamma": self.gamma}
221        if not return_surrogates:
222            results.pop("surrogates")
223        return results
224    def mcssa_procrustes(self,up_perc,down_perc,ns=100000,max_iter=1000,return_surrogates=False):
225        """
226        Procrustes-aligned MCSSA.
227
228        For each surrogate:
229        - compute surrogate covariance eigenvectors and eigenvalues,
230        - scale surrogate EOFs by sqrt(eigenvalues),
231        - compute the orthogonal Procrustes transform that best aligns the
232          surrogate-scaled EOFs with the data-scaled EOFs,
233        - apply the transform to surrogate eigenvalues and collect projected
234          diagonal spreads to form percentile confidence bounds.
235
236        Args:
237            up_perc (float): Upper percentile for confidence.
238            down_perc (float): Lower percentile for confidence.
239            ns (int, optional): Number of surrogate realizations.
240            max_iter (int, optional): Max iterations for AR(1) estimation.
241            return_surrogates (bool, optional): If True include surrogates in the result.
242
243        Returns:
244            dict: {
245                "data_eigenvalues": ndarray,
246                "upper_confidence": ndarray,
247                "lower_confidence": ndarray,
248                "surrogate_eigenvalues": list of ndarrays,
249                "surrogate_eigenvectors": list of ndarrays,
250                "procrustes_transformations": list of ndarrays,
251                "surrogates": ndarray,
252                "alpha": float,
253                "gamma": float
254            }
255
256        Raises:
257            ValueError: If ns < 1 or invalid percentile arguments.
258        """
259        if ns < 1:
260            raise ValueError("ns must be >= 1")
261        if not (0 <= down_perc <= 100 and 0 <= up_perc <= 100):
262            raise ValueError("Percentiles must be between 0 and 100")
263        if down_perc >= up_perc:
264            raise ValueError("down_perc must be less than up_perc")
265        self.gamma=round(self.gambar(max_iter=max_iter),8)
266        self.alpha=round(self.alph(max_iter=max_iter),8)
267        surrs=[]
268        for i in range(ns):
269            surrs.append(self.ar1_model(self.alpha,self.gamma))
270        surrs=np.vstack(surrs)
271        surrs=surrs-np.mean(surrs,axis=1,keepdims=True)
272        self.surrs=surrs
273        l_rs, e_rs, e_rscaled=[],[],[]
274        for i in range(np.shape(surrs)[0]):
275            T_r=Tmat(surrs[i,:],self.m)
276            c_r=(1/self.nd)*np.transpose(T_r)@T_r
277            l_r, e_r=np.linalg.eigh(c_r)
278            idx=np.argsort(l_r)[::-1]
279            l_r=l_r[idx]
280            e_r=e_r[:,idx]
281            sigma_r=np.sqrt(l_r)
282            e_scaled=e_r*sigma_r[np.newaxis,:]
283            l_rs.append(l_r)
284            e_rs.append(e_r)
285            e_rscaled.append(e_scaled)
286        T_dat=Tmat(self.data-np.mean(self.data),self.m)
287        c=(1/self.nd)*np.transpose(T_dat)@T_dat
288        l,e=np.linalg.eigh(c)
289        idx=np.argsort(l)[::-1]
290        l=l[idx]
291        e=e[:,idx]
292        sigma=np.sqrt(l)
293        e_scaled=e*sigma[np.newaxis,:]
294        T_es=[]
295        for e_r_sig in e_rscaled:
296            M=e_r_sig.T@e_scaled
297            U,S,Vt=np.linalg.svd(M)
298            T_e=U@Vt
299            T_es.append(T_e)
300        l_proj_list=[]
301        for l_r,T_e in zip(l_rs,T_es):
302            l_rmat=np.diag(l_r)
303            l_proj=T_e.T@l_rmat@T_e
304            l_proj_diag=np.diag(l_proj)
305            l_proj_list.append(l_proj_diag)
306        l_projs=np.vstack(l_proj_list)
307        ups=np.percentile(l_projs,up_perc,axis=0)
308        downs=np.percentile(l_projs,down_perc,axis=0)
309        results={"data_eigenvalues": l,
310                   "upper_confidence": ups,
311                   "lower_confidence": downs,
312                   "surrogate_eigenvalues": l_rs,
313                   "surrogate_eigenvectors": e_rs,
314                   "procrustes_transformations": T_es,
315                   "surrogates": surrs,
316                   "alpha": self.alpha,
317                   "gamma": self.gamma}
318        if not return_surrogates:
319            results.pop("surrogates")
320        return results